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I.  INTR0DLCT10N 


The  analysis  and  synthesis  of  linear  feedback  control 
systems,  or  the  compensation  of  same,  can  be  realized  by 
three  general  methods.  The  first  of  these  can  be  called  the 
integral  method.  Given  a  control  system,  described  by  a  set 
of  differential  equations,  one  selects  a  cost  function  to  be 
minimized  with  respect  to  certain  variable  system  parameters. 
The  major  drawback  with  this  method  is  the  difficulty  of 
varying  more  than  one  parameter  at  a  time.  The  second 
method  is  the  Bode  frequency  response  technique  whereby  the 
system's  open  loop  transfer  function  is  manipulated  to  obtain 
the  desired  system  response.  This  method  also  has  its 
inherent  weaknesses:  difficulty  of  application  to  non-unity 
feedback  control  systems,  difficulty  in  interpreting  the 
closed  loop  transient  response  in  terms  of  the  open  loop 
frequency  response,  and  difficulty  of  varying  more  than  one 
parameter.  Third  are  the  algebraic  methods.  Within  this 
category  can  be  included  the  familiar  root  locus  method. 

Here,  a  graphical  technique  is  provided  by  which  the  set  of 
all  points  which  could  potentially  be  made  roots  are  plotted 
in  the  S-plane.  The  root  locus  method  is  a  valuable  and 
powerful  tool  when  only  one  parameter  is  varied;  results 
are  less  satisfactory  for  two  parameters  and  of  little  use 
when  three  or  more  parameters  are  involved. 


Methods  for  studying  the  parameter-root  relationship 
when  two  or  more  parameters  are  variable  are  clearly  of 
considerable  value.  For  a  linear  system,  the  set  of 
differential  equations  that  describe  that  system  can  be 
transformed  into  algebraic  equations  and  manipulated  to 
provide  a  characteristic  polynomial.  Since  the  coefficients 
of  the  characteristic  polynomial  are  deterimined  by  the 
system  parameters,  it  follows  that  some  relationship  exists 
between  the  value  of  any  parameter  and  the  value  of  the 
characteristic  roots.  In  reference  (1),  Mitrovic  developed 
an  algebraic/graphical  method  for  obtaining  the  roots  of  a 
polynomial  in  terms  of  two  variable  parameters.  In 
references  (2),  (3),  and  (4),  Choe ,  Hyon ,  and  Nutting, 
respectively  developed  and  extended  the  Mitrovic  method  to 
the  compensation  of  linear  continuous  feedback  control 
systems.  The  disadvantage  of  the  Mitrovic  method  is  that  the 
variable  parameters  may  appear  in  no  more  than  two  coeffi¬ 
cients  of  the  characteristic  equation,  which  limits  the 
flexibility  of  the  technique.  In  reference  (5),  Siljak 
introduced  a  method  for  obtaining  the  roots  of  a  polynomial 
in  terms  of  two  variable  parameters  that  may  appear  in  any 
and  all  of  the  coefficients  of  the  polynomial.  Later,  Thaler 
and  Towill  [Ref.  6]  extended  this  method  to  the  compensation 
of  linear  continuous  feedback  control  systems.  It  is  from 
the  latter  work  that  the  ensuing  parameter  plane  equations 


were  developed.  General  methods  of  compensation  will  be 
presented,  and  an  attempt  will  be  made  to  relate  the  root 
locus  and  the  parameter  plane  methods  as  a  set  of  comple¬ 
mentary  techniques  which,  when  applied  in  tandem,  represent 
the  most  satisfactory  tool  to  date  for  designing  linear 
feedback  control  systems. 

The  parameter  plane  method,  which  works  well  for  two 
variable  parameters  and  which  may  be  extended  to  three  or 
more  parameters,  is  purely  algebraic,  and  the  resulting  plots 
are  valuable  aids  to  analysis.  The  term  parameter  plane 
comes  from  the  plot  for  two  parameters--in  a  rectangular 
coordinate  space  one  parameter  will  define  the  abscissa  while 
the  second  parameter  defines  the  ordinate  (the  S-plane  is 
inconvenient  for  presenting  the  desired  results).  Three 
parameters  define  a  3-dimensional  parameter  space,  etc.  For 
design  problems  it  is  convenient  to  think  of  the  algebraic 
calculations  as  a  mapping  procedure.  By  choosing  a  point  on 
the  S-plane,  the  characteristic  polynomial  acts  as  a  mapping 
function  whereby  the  point  may  be  "mapped"  onto  the  alpha- 
beta  plane  (alpha  and  beta  are  the  two  variable  parameters 
to  be  used  throughout  the  remainder  of  this  text).  The 
relationship  between  being  able  to  place  the  roots  of  a 
polynomial  at  specific  locations  in  the  S-plane  and  the 
compensation  of  linear  feedback  control  systems  is  as 
follows.  A  feedback  control  system,  including  any  added 
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compensators  which  may  contain  variables,  can  be  reduced  to  a 
ratio  of  two  polynomials  (the  closed  loop  transfer  function). 
A  specified  system  response  in  terms  of  overshoot,  bandwidth, 
settling  time,  steady-state  accuracy,  etc.,  can  theoretically 
be  obtained  by  placing  a  pair  of  complex  conjugate  roots  of 
the  characteristic  equation  at  a  specific  location  in  the  S- 
plane,  while  ensuring  that  the  real  part  of  this  complex  root 
pair  (the  dominant  roots)  is  smaller  in  magnitude  than  the 
real  parts  of  the  remaining  roots  of  the  characteristic 
equation.  The  problem  of  compensation,  thus  of  feedback 
control  system  design,  reduces  to  one  of  placing  the  dominant 
roots  of  the  characteristic  equation  at  the  desired  location. 
The  ability  of  the  parameter  plane  method  to  achieve  this 
goal  will  become  obvious. 
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II.  DERIVATION  OF  PARAMETER  PLANE  EQUATIONS 


A  linear  feedback  control  system's  characteristic  equation 
can  be  expressed  as  a  polynomial  of  the  following  form: 


f(s)  =  2.  a,S  =  0,  where 

k=0  K 


(2-1) 


a^  ( k=0 , 1 , . . . ,m)  are  real  coefficients 

2 

S  =  —  (T +j u)  =  -|oj+  jw  /l-£ 

is  the  undamped  natural  frequency  and 
is  the  relative  damping  coefficient 
In  reference  (5)  it  is  noted  that  S  may  be  represented  by 
the  following: 


Sk  =  u;k(Tk(-£)+j  /  1-|2  uk(-o) 


(2-2) 


where 

Tk(-C)  =  (-l)kTk(U  and  uk(  -  «)“( -1  )k+1UR(  |)  . 

T  (£)  and  U  (£)  are  Chebishev  functions  of  the  first  and 
K  K 

second  kind  respectively.  Values  of  zeta  and  omega  will  be 
considered  such  that  0<|<1  and  0<u><oo.  Values  of  Tk  and  Uk 
are  tabulated  in  .various  appendixes.  More  important  to  digital 
computer  analysis,  they  can  be  obtained  from  the  following 
recursive  relations: 


9 


Tk+i(U  -  2V*}  +  Tk-i(*}  =  0 


uk+l(?)  -  2uk(^  +  uk-l(?)  =  0 


(2-3 


Here,  T0(|)=l,  T^f)-?,  UQ(n  =  0,  U1(g)  =  l.  Substituting 
equation  (2-2)  into  (2-1)  and  setting  the  real  and  imaginary 
parts  to  zero  independently,  one  obtains: 


S  a.u*T.(-5)«0 
k=0  K  K 


Z  a  cj^U  (-0=0 
k=0  K  K 


(2-4 


Employing  equations  (2-3),  one  obtains  from  equation  (2-4): 


Z  (-l)ka,cAj,  (*)  =  0 
k=0  K 


k  k 

z  (  -1  )Ka,  wU  (  5  )  =  0 
k=0  K  K 


(2-5 


Now  consider  the  coefficients  a,  of  the  characteristic 

k 

equation  (2-1)  as  linear  functions  of  the  variable  system 
parameters,  a  and  (3 ,  as  follows: 


a.  =  b,  a  +  c,  B  +  d, 


(2-6 


L’sing  this  relation  for  a^,  equations  ( 2-5 )  become: 


a B  +  SC1  +  D1  =  0 


(2-7) 


*  6C2  +  D2  "  0 


where 
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l  <-i>k  tvk  uk_1 


B, 


m  k  k 

£•  (-1)K  b,  <3  U 

k=0  K 


„  ”  ,  1Ak  k 

C  =  2  (-D  c.w 

1  k=0 


k-1 


m  ]•  k 

s  <-1)  V  \ 


(2-8) 
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m 

Z 

k=0 


(-1)k  ^  Uk- 


°2  = 


m 

r 

k=0 


,  ,k  ,  k 
(-1)  dkw 


Since  equations ( 2-7 )  are  linear  in  the  two  unknowns  alpha  and 
beta.  Cramer's  rule  may  be  applied  to  obtain: 


Equations (2-9)  are  now  functions  of  zeta  and  omega.  Hence, 
by  fixing  either  zeta  or  omega  and  varying  the  remaining 
parameter,  the  constant  omega  or  constant  zeta  S-plane 
contours  respectively  can  be  mapped  into  the  real  domain 
of  the  alpha-beta  or  parameter  plane. 

In  reference  (5)  the  following  relationships  are  noted: 

Sk  =  Pk  +  jo>/  I-*2  Qk 

P  +  2ojP  +  o)Jp  =0 
k+1  k  k-1 

Qk+1  +  2o;Qk  +  w2Qk_1  =  0  (2-10) 

Pk  >  -crtQk  -  u,2Qk_1 


Po  =  P1  '  -«">  «1  '  ! 

Pk  and  Qk  are  related  to  the  Chebishev  functions  by: 

Pk  =  wkTk(-f)  =  (-l)ku'kTk(5) 

=  Wk~1Uv(-4)  =  ( -1 )k+1wk-1u. (£  ) 


(2-11) 


Q, 


By  employing  equations  (2-10)  and  (2-11),  one  obtains 
(proceeding  as  before); 
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k=0 


ak^k-l 
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k=0 


ak*^k 


0 


(2-12) 


Combining  equations  (2-6)  and  (2-12)  with  Cramer's  rule,  one 
again  arrives  at  equations  ( 2-9 )  where  the  following  expressions 
now  apply: 
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(2-13) 
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Equations  (2-9)  and  (2-13)  are  useful  for  mapping  constant 
zeta-omega  curves  from  the  S-plane  into  the  parameter  plane. 

As  will  be  demonstrated  later,  these  curves  play  an  important 
role  in  dominance  considerations. 

If  the  complex  variable  S  is  substituted  in  equation  (2-1) 
by  letting  S  =— a,  where  sigma  corresponds  to  values  of  S  along 


v; 

y.< 


%  , 
% 


the  real  axis,  then  according  to  equation  (2-6)  the 
characteristic  equation  (2-1)  becomes: 

m  k  k  m  kkm  kk 

a  2  (-1)  b.a  +  0  2  (-1)  c .a  +  L  (-1)  ci  crK  =  0 

k=0  K  k=0  K  k=0  K 

(2-14) 

The  above  expression  represents  a  straight  line  in  the  alpha- 

beta  plane  for  a  given  value  of  sigma.  Hence  a  point  on  the 

real  axis  in  the  S-plane  maps  into  a  straight  line  in  the 

alpha-beta  plane.  In  addition,  for  given  values  of  alpha, 

beta,  and  sigma  which  satisfy  equation  (2-14),  the 

characteristic  equation  (2-1)  must  have  a  real  root  at  minus 

sigma.  On  the  constant  zeta  and  omega  curves  previously 

defined,  for  certain  values  of  alpha  and  beta  (say,  for 

values  obtained  from  equations  (2-9)  for  given  values  of 

zeta  and  omega)  the  characteristic  equation  will  have  a  pair 

— . 

of  complex  conjugate  roots  at  S  =  -§o»  +  j w/  1-C 

The  significance  of  the  above  discussion  is  that  by 

applying  equations  (2-9)  and  (2-14)  one  can,  for  a  specified 

value  of  zeta,  omega,  and  sigma,  compute  the  value  of  alpha 

and  beta  such  that  the  characteristic  equation  will  have  a 

2 

pair  of  roots  at  S  =  *  The 

remaining  roots  of  the  characteristic  equation  can  then  be 
calculated  by  dividing  out  the  two  known  or  specified 
roots.  This  method,  where  zeta,  omega  and  sigma,  or  simply 
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zeta  and  omega  are  specified,  and  where  the  computations  for 
alpha  and  beta  are  done  algebraically,  will  be  referred  to  as 
the  algebraic  parameter  plane  solution. 

To  solve  the  problem  in  general  for  all  values  of  zeta, 

omega,  and  sigma,  it  becomes  necessary  to  plot  a  family  of 

parameter  plane  curves  for  various  values  of  zeta,  omega, 

sigma,  and  if  desired,  zeta-omega.  On  the  resulting 

parameter  plane  plot  one  can,  by  choosing  an  operating  point, 

graphically  read  from  the  curves  the  values  of  alpha  and  beta 

and,  hence,  the  values  of  the  nr.  roots  of  the  corresponding 
t  h 

m  order  characteristic  equation.  This  latter  method  will 
be  referred  to  as  the  graphical  parameter  plane  solution. 

The  algebraic  solution  has  the  advantage  that  the  labor 
of  plotting  the  curves  can  be  avoided,  but  the  disadvantage 
remains  that  without  the  curves  it  is  difficult  to  pick  the 
optimum  values  of  zeta  and  omega  so  as  to  ensure  dominance 
while  still  meeting  the  system  specifications.  The  graphical 
solution  has  the  advantage  that  one  has  a  "picture"  of  the 
way  the  characteristic  roots  move  about  in  the  S-plane  as 
alpha  and  beta  are  varied.  This  enables  one  to  choose  the 
values  of  alpha  and  beta  corresponding  to  the  best  values  of 
zeta,  omega,  sigma,  and  zeta-omega  for  all  roots  of  the 
characteristic  equation.  This  feature  of  the  parameter 
plane  points  out  a  strong  justification  for  attempting  to 
obtain  the  parameter  plane  curves.  And  with  the  employment. 


of  a  digital  computer  and  an  appropriate  algorithm  to  realize 
the  parameter  plane  curves,  the  advantage  of  the  algebraic 
method  becomes  muted. 

Recursion  methods  (equation  (2-3))  are  by  no  means  the 
only  methods  of  producing  algebraic  and  graphical  parameter 
plane  data.  Thaler  and  Karmarkar  [Ref.  7]  describe  a  matrix 
solution  to  the  parameter  plane  problem.  Essentially,  a 
matrix  of  coefficients  may  be  manipulated  to  obtain  the 
following  general  form: 


1(5 


'  vs* 

.*  «• 
■  »  ' 


V.  rrv 

cs;- 


V 


■  W  ja  ^  f  a  .  .Jl 


where  b^,  c^,  d^,  and  are  as  described  before,  and: 


e  (k=l,  m)  are  the  coefficients  associated  with  the 

non-linear  alpha-beta  product  terms 

__  ° 

R..  is  the  sum  of  the  m-2  roots  of  the  polynomial  character¬ 
istic  equation  taken  one  at  a  time 


m-2 

R  „  is  the  sum  of  the  m-2  roots  taken  m-2  at  a  time 
m-2 


Further,  by  the  application  of  appropriate  row  operations, 
this  matrix  may  be  reduced  to  the  following  row-echelon  form: 


1 

0 

0 

• 

• 

0 

kll 

k12 

a 

0 

1 

0 

• 

• 

0 

k2 1 

** 

to 

to 

6 

0 

0 

1 

• 

• 

0 

k31 

k32 

m-2 

Rm  9 
m-2 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

1 

0 

0 

0 

• 

• 

1 

k  , 
ml 

km2 

aS 

=  0 


(2-15) 
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For  the  case  when  all  coefficients  of  the  characteristic 
equation  are  linear,  i.e.  e^(k=l , .  .  .  ,ra)  =  0,  then 
( k=l , . . . ,m )  =  0 ,  and 


a 


6  =  -K 


21 


as  obtained  from  the  first  two  rows  of  the  matrix  equation. 

One  should  note  that  in  arriving  at  equations  (2-15)  , 

2 

approximately  m  row  operations  are  required  for  the  row- 
echelon  matrix  formulation  for  each  point  of  the  parameter 
plane  curves  (e.g.,  each  time  either  zeta  or  sigma  are 
varied).  Compare  this  with  the  approximate  m  calculations 
required  to  obtain  the  recursion  equations  of  the  previous 
chapter,  and  the  matrix  method  becomes  relatively  inefficient 
for  larger  order  systems. 

One  should  not,  however,  discard  the  matrix  approach 
entirely.  For  small  order  characteristic  equations,  this 
technique  compares  favorably  with  the  recursion  method. 

And  when  the  variable  parameters  are  non-linear — when  one 
must  deal  with  alpha-beta  product  terms--the  matrix  approach 
affords  a  more  direct  method  of  obtaining  the  alpha-beta 
pairs.'*’  Whether  the  recursion  or  matrix  method  is  utilized 


From  equations  (2-15),  one  obtains  the  two  quadratic  forms 

K22°2-(K21  K12-K11.K22-1,“+K11-° 

K126  *<K21  K12-KU  K22tl)8+K21=° 
from  which  alpha  and  beta  are  easily  derived. 


•«<- 


18 


'  -.'W 


III.  APPLICATION  OF  THE  PARAMETER  PLANE  METHOD 

A.  ALGEBRAIC  SOLUTION 

In  this  section  it  will  be  assumed  that  the  system 
performance  specifications  can  be  met  by  placing  a  pair  of 
complex  conjugate  roots  at  a  specific  location  (i.e.,  by 
choosing  appropriate  values  for  zeta  and  omega).  If  after 
computation  of  the  necessary  values  of  alpha  and  beta  to 
locate  the  roots  as  desired  it  is  found  that  these  specified 
roots  are  not  dominant,  then  either  a  different  value  of  zeta 
and/or  omega  must  be  used  (possibly  at  the  sacrifice  of  some 
performance  measure),  or  a  different  method  of  compensation 
will  have  to  be  attempted.  In  a  later  section  a  method  will 
be  addressed  whereby  the  dominancy  requirement  may  be 
achieved. 

1 .  Feedback  Compensation 

For  a  unity  feedback  control  system,  let 


G 


K 

e(S) 


K 


am  qIti  1 ,  ,  ql 

S  +e  „  b  + .  .  .  +eT  b 
m-1  L 


(3-1 


where  K  is  the  forward  path  gain  (a  variable)  and  e(S)  is  a 
polynomial  in  S  representing  the  poles  of  the  open  loop 
transfer  function  of  the  uncompensated  system.  In  equation 


Si 


(3-1),  L  corresponds  to  the  system  type--for  a  type  0 
system,  L=0 ,  for  type  1,  L=l,  etc.  The  system's  error 
coefficient  is  defined  as: 

K  =  lim  SLG  (3-2) 

e  S-0 

where  G is  t'he  open  loop  transfer  function  of  the 
compensated  system. 

a.  Tachometer  Plus  Acceleration  Feedback 

In  order  to  achieve  the  system  performance 
specifications,  a  feedback  compensator  must  be  introduced. 
Let 

H  =  K, S  +  K  S2 
t  a 

The  resulting  compensated  system's  characteristic  equation 
becomes: 

e(S)  +  K( K  S+K  S2)  =  0  (3-3) 

L  3. 

and  by  expanding  e(S),  equation  (3-3)  becomes: 

Sm+e  1Sm-1+. . .+(e  +KK  )S2+( e +KK )S+e  +k  =  0  (3-4) 

m-1  2  a  1  t  o 

where  L  is  zero  for  a  type  0  system  (the  most  general  case). 
The  following  result,  also  apply  to  a  type  1  system  if  eQ  is 
set  to  zero,  and,  similarly,  for  a  type  2  system  if  both 
eQ  and  e^  are  set  to  zero,  etc.  Combining  equations  (3-2), 
(3-3),  and  (3-4)  the  error  coefficient  becomes: 


•  *'  »  4  -Jl  >  »  »  .■  •  -J*  'J  J  V  "VJ*  «T  *"  .T 


e  S-0  e(  S  )+K(  K  S+K  S*)  co  '  “ 

L  3. 

or  for  a  type  1  uncompensated  system: 

K  =  K 

Ke  -  e1+KKt 

or  for  a  type  2  uncompensated  system: 

K  -  K 

e  "  KKt 


(3-6 


(3-7 


Note  that  if  the  uncompensated  system  is  type  2 ,  the 
compensated  system  would  be  type  1  if  tachometer  feedback 
or  tachometer  plus  acceleration  feedback  is  used. 

In  the  compensated  system's  characteristic 
equation  (3-4)  let  alpha  =  KK  and  beta  =  KK  .  Equation 

cl  L 

(3-4)  then  becomes: 

Sm  +  e  .  Sm_1  +. . .  +  (e0  +  a)S2  +  (e  +6)S+e  +K  =  0 
m-1  3  1  o 

Recalling  equation  (2-6)  where  in  general  the  coefficients 
f  the  characteristic  equation  are  of  the  form: 


V+ck6+dk 


and  letting  m=k,  then  from  equations  (2-8)  one  obtains: 


B1  =  (-!)-»'  U  = 


B2  =  *  U2 


c  =  -«u  =  0 

1  o 


c2  =  -«Uj  =  -• 


(3-8) 


m  .  ,  m  , 

D  =  Z  (-l)Kd  .<*>KU  D  =  z  ( -1  )Kd  >“  kU 

1  k=0  k  k'1  2  k=0  k  k 


since  Uq=0  and  U^  =  l.  From  equations  ( 2-9 )  one  derives: 

m  k 

C  D_-C9D1  (_1)  dkUk-l  m  .  ,  0 

a  =  rV-B-C  =  1 -  --X  (  -1  )kd.  «>k_2U 

B1C2-B2C1  -c?  k=0  k 


k-1 


(3-9) 


m 


8  '  k(o  <-1)Vk'1<VU20k-l> 


If  alpha  and  beta  are  linear  functions  of  K,  the  forward 

path  gain,  one  can  use  the  steady-state  error  specification 

to  define  K  in  terms  of  alpha  and/or  beta.  Since  zeta  and 

omega  were  assumed  to  be  specified,  then  from  equations 

(3-9)  one  may  solve  for  alpha  and  beta.  From  this,  K  and 

o 

K  are  readily  determined, 
t  J 

Example  3-1 


The  system  of  Figure  (3-1)  is  to  be  compensated  by  using 
tachometer  plus  acceleration  feedback .  The  system 
specifications  are  as  follows: 


Figure  3-1 


From  which  K_>12+6KK^ .  The  compensated  characteristic 
equation  is 


S3+S2(3+KK  )+S( 2+KK  )+K  =  0 


(3-10) 


Letting  alpha  =  KK  and  beta  =  KK  equation  (3-10)  becomes: 

3  "t 

S3+-S2(3+cv)+S(2+^)+K  =  0  (3-10a) 

From  equations  (3-8): 

B  =  100 
C1  =  ° 

D  =  -1100-K 

and  from  equations  (3-9): 

„  -  10( -1100-K )  140(-1100-K)+56000  ,, 

-1000  ’  P  -1000 

From  the  steady-state  accuracy  specifications,  then,  it  is 

necessary  that  K>12+6/?;  let  K=12+6/?.  From  equations  (3-11) 

it  is  found  that  8=623,  hence  K=3750.  Therefore,  a=48.5, 

and  since  <*=KK  and  #=KK.  : 

a  t 

48  5 

K  =  ~ ^  =  0.0129 
a  3750 


B2  =  140 

C2  =  _1° 

D  =  -1120 


K, 


623 


0. 1661 


The  compensated  system's  characteristic  equation  becomes 


S3+51 . 5S2+625S+ 3750=0  (3-12) 

2 

Now  zeta=0.7  and  oo  =  10  corresponds  to  S  +14S+100=0. 

Dividing  equation  (3-12)  by  this  quadratic,  the  remainder 
is  S+37.5.  Since  zeta*omega  of  the  desired  roots  =  7<<37.5, 
the  complex  roots  are  dominant  and  the  problem  is  solved, 
b.  Tachometer  Feedback  Only 

Let  H  =  K^S.  The  characteristic  equation  of  the 
compensated  system  becomes: 

Sm+em_iSin_1  +  ..  .+e2s2+(ei+a)s+e0+3=° 

Proceeding  as  in  the  previous  example,  one  obtains: 


Bi  -  0  B2  ~ 

C1  =  -1  c2  =  0 

m  k  k  m  k  k 

D,  =  Z  (-l)V  /U,_i  D_  =  r  (-l)Kd^KUv 
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For  a  specified  value  of  zeta  and  omega,  alpha  and  beta  can 
be  obtained  from  equations  (3-13).  The  error  coefficient  is 
then  determined  directly  from  equations  (3-5),  (3-6),  or 
(3-7).  Thus  the  error  coefficient  is  fixed  for  a  given 
value  of  zeta  and  omega,  and  if  this  parameter  is  to  be  met, 
the  values  of  Zeta  and  omega  may  require  adjustment.  One 
possible  approach  might  be  to  fix  zeta  at  some  value, 
whereby  from  the  given  and  equations  (3-5),  (3-6),  or 
(3-7)  alpha  could  be  computed.  Equations  (3-13)  could  then 
be  solved  for,  first,  omega  and  then  beta.  The  calculations 
would  prove  tedious,  however. 


The  same  system  as  used  in  example  (3-1)  will  be  studied 
here,  this  time  with  tachometer  feedback  alone.  The  same 
system  performance  specifications  are  to  be  met,  namely, 

Kg>6 ,  zeta  =  0.7,  and  omega  =  10.  The  compensated  system's 
characteristic  equation  becomes: 

S3+3S2+(2+KKt)S+K  =  0  (3-14) 

Letting  a^KK^  and  j3=  K  here,  equation  (3-14)  becomes: 

S3+3S2(2  +  or)S+/3  =  0 

From  equation  (3-13)  it  is  found  that: 

a  =  -2+30(1. 4)-100(0. 96)  =  -56 

Since  alpha  is  negative,  it  is  seen  that  positive  tachometer 

feedback  is  required.  Further,  it  is  found  that  the 

remaining  root  (when  equation  (3-14)  is  divided  by 
2 

S  +14S+100)  is  positive;  the  system  is  unstable.  Hence 

the  desired  system  specifications  cannot  be  met  with 

tachometer  feedback  alone. 

c.  Accelerat ion  Feedback  Only 
2 

Let  H  =  K^S  .  The  characteristic  equation  of 
the  compensated  system  becomes: 


Proceeding  as  before,  where  now  a=KK  and  /3= K: 

a. 


Bi  =  u,2u1=u’2 


B2  =  •“  U2 


C,  = 


U-i  =  -1 


C„  =  0 


m 


D.  = 


k=0  K  K  x 


m 


= 


r  (-i)kd  -ku 
k=0  K  k 


Solving  for  alpha  and  beta  yields: 


“Dr 


a  ~ 


2n 
u  Ur 


_  i_  (_l)kd  u>  k_2U 
u2  ^  ak  Uk 


(3-15) 


a  v  r  \k ,  kTT  1  ™  .  .k  ,  kTT 

1  (-1)  dk-  U  "  ip  1  (-1)  dku  Uk 

k=0  1  2  k=0  K  K 


Calculations  for  alpha,  beta,  and  Kg  are  performed  in  the 
same  manner  as  with  the  preceding  tachometer  feedback 
example . 

Example  3-3 

The  same  system  of  examples  (3-1)  and  (3-2)  will  now  be 
compensated  using  acceleration  feedback  alone.  As  before, 

K  >6,  zeta=0.7,  and  omega=10.  Therefore  K  =—  and  the  error 


coefficient  is  unaffected  by  the  acceleration  feedback.  Hence 
one  can  conveniently  choose  K=12  to  meet  the  specifications. 
The  compensated  system's  characteristic  equation  becomes: 

S3+( 3+KK  )S2+2S+K  =  0  (3-16) 

a 

If  K  in  equation  (3-16)  is  set  equal  to  12  as  prescribed,  only 
one  parameter  remains  and  the  parameter  plane  equations 
produce  an  indeterminate  solution.  If  K  is  left  as  the 
variable  beta,  then  equation  (3-16)  becomes  (after  the  usual 
subst itut ions  ) : 

S3+(3+a)S2+2S+/3  =  0 

By  employing  equations  (3-15),  one  obtains  a=4  and  /3=  -700. 
Since  beta  is  negative,  it  is  concluded  that  the  desired 
roots  (i.e.,  desired  values  of  zeta  and  omega)  cannot  be 
realized  using  acceleration  feedback  alone,  and  of  course 
neither  can  the  desired  error  specification  be  obtained. 

One  would  therefore  choose  an  alternate  method  of 
compensat ion . 

If  one  chooses  to  use  feedback  compensation  then 
perhaps  tachometer  plus  acceleration  feedback  might  be 
attempted  first  using  equations  ( 3-9 )  and  the  appropriate 
steady-state  error  specification.  If  the  specif ications  cannot 
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be  met  in  this  manner,  then  it  follows  that  neither 


t  acLometer  nor  acceleration  fL-cdbmk  alone  will  Mill  nc,  In 
this*  case  either  the  system's  specifications  must  be  eased  or 
another  type  of  compensation  must  be  utilised.  If  it  is 
found  that  the  specifications  are  achievable  with  the  combined 
tachometer  and  acceleration  feedback,  then,  if  desired, 
equations  (3-13)  and  (3-15)  can  be  employed  t  *  j  investigate 
the  feasibility  of  tachometer  or  acceleration  feedback  alone, 
respect i ve  ly . 

d.  Case  For  Which  Feedback  Is  Not  Available  Near 
The  Forward  Path  Amplifier 


Figure  3-3 

Figure  3-3  shows  a  system  similar  to  that  used 
in  example  (3-1)  except,  that  now  the  feedback  is  inserted  at 
the  output  terminals  of  the  amplifier  represented  by  gain  K. 
This  illustrates  a  system  for  which  it  may  not  be  possible 
or  practical  to  access  the  input  terminals  of  the  error 
detector.  This  problem  will  be  solved  by  means  of  an  example. 


Example  3-4 

As  before,  the  same  system  specifications  are  to  be  met, 
i.e.,  K^>6,  zeta=0.7,  and  omega=10.  The  characteristic 
equation  becomes: 

S3+( 3+K  )S2+( 2+K, )S+K  =  0 
a  t 

Letting  or=K  and  (i=K  : 

S3+(3+a)S2  +  (2  +  |3)S+K  =  0  (3-17) 

Comparison  of  equations  (3-17)  and  (3-10a)  show  that  they 
are  identical,  that  is,  the  solution  obtained  for  alpha  and 
beta  in  example  (3-1)  applies.  There,  alpha  and  beta  were 
found  to  be  48.5  and  623,  respectively,  while  K=3750.  For 
the  present  example  no  further  computations  are  necessary 
to  find  K  and  K  ,  since  they  are  now  the  parameters  alpha 
and  beta.  This  points  out  an  important  advantage  of  the 
parameter  plane  method,  namely,  that  the  solutions  depend  only 
on  the  characteristic  equation  and  not  on  the  system  from 
which  the  characteristic  equation  was  formed.  This  principle 
can  similarly  be  applied  to  control  problems  involving 
tachometer  or  acceleration  feedback  alone. 

2 .  Cascade  Compensation 

For  a  unity  feedback  control  system  let  G  have  the 


form  of  equation  (3-1): 


where  again  K  is  the  forward  path  gain  (a  variable)  and 
e(s)  is  a  polynomial  in  S  representing  the  poles  of  the 
open  loop  transfer  function  of  the  uncompensated  system. 
The  letter  L  again  indicates  the  system  type.  If  in  order 
to  satisfy  the  system's  requirements  a  cascade  compensator 
is  required,  then  let: 

r  =  P(S+Z) 

C  Z(S+P) 


With  a  d.c.  gain  of  unity,  placement  of  this  compensator  in 
the  forward  path  will  not  affect  steady-state  accuracy. 

With  Gc  as  indicated  here,  the  values  of  P  and  Z  are 
computed  to  obtain  the  desired  system  response.  If  P  is 


less  than  Z,  a  lag  network  is  required  and  the  factor  of 
P 


of  the  compensator  is  inherently  present  due  to  the 
physical  nature  of  the  compensator  (usually  an  R-C  network) 
In  this  case  all  forward  path  amplifier  gains  can  remain 
unaltered  to  meet  the  specified  accuracy  demands.  If, 
however,  the  computed  value  of  P  is  greater  than  Z,  a 


lead  network  is  required  and  the  compensated  system's 


forward  path  gain  must  be  raised  by  the  factor  of  ^ 


meet  the  accuracy  specifications.  As  the  physical  nature 


P  . 


of  the  lead  network  is  such  that  the  factor  77  is  not 


> .-v 


inherently  present,  this  factor  must  be  provided  either  by 
adding  an  amplifier  in  cascade  with  the  lead  network  or  by 
raising  the  gain  K  of  the  existing  amplifier  as  required  to 
achieve  steady-state  accuracy. 

Continuing,  the  compensated  system's  forward  path 
transfer  function  is: 

r  _  r  r  _  K  P  S+Z  _  V(S+v)  K 

CC  C  e(s)  '  Z  *  S+P  S+P  ’  e(s) 


Applying  the  definition  of  the  error  coefficient  one 
obtains : 


K 

e 


lim  S^‘ 
s  -*-0 


K 


e(s) 


Y(s+y) 

(S+P) 


K 


e 


L 


and  again  assuming  a  type  0  system  where  L=0,  the  compensated 
system's  characteristic  equation  becomes: 


e(s)(S+P)+KY(S+y)  =  0 


or  after  expansion: 

Sra+1+(p+e  )Sm+( Pe  +e  „)Sm_1+... 
m-1  m-1  m-2 

+(Pe„+e. )S2+(KY+Pe  +e  )S+P(e  +K)  = 


0 


(3-18) 


Letting  =P  and  3  =  y,  equation  (3-18)  becomes: 


Sm+1  +  (a+e  )Sra+(e  <*  +  em  9)Sm  X+... 

m-l  m-1  m-2 

+  (e2u+e1)S2  +  (K3+e1ot+e0)S+a(e0+K)  =  0  (3-19) 


Comparison  of  equation  (3-19)  with  the  general  form  of  the 
characteristic  equation  as  specified  in  equat ions ( 2-1 )  and 
(2-6),  it  is  apparent  that  K=m+1 ,  bQ=e0+K,  c0=dQ=0,  b1=e1 , 
c  ^ =K ,  d-^=eQ,  b2=e2  ,  c2  =  0,  d2=e^,  etc. 

It  is  important  to  note  that  the  parameter  plane 
variable  beta  represents  the  pole-to-zero  ratio  of  the 
cascade  compensator.  The  S-plane  can  be  divided  into  regions 
where  lag  compensation  or  lead  compensation  is  needed.  By 
mapping  of  variables  in  the  above  manner,  the  parameter 
plane  can  effectively  be  divided  into  corresponding  regions 
above  and  below  the  straight  line  3=1.  Then,  for  values  of 
beta  less  than  one  a  lag  network  is  required  and  for  beta 
greater  than  one  a  lead  network  is  needed.  In  addition,  if 
beta  is  less  than  0.1  or  greater  than  10,  a  multiple  lag  or 
multiple  lead  network,  respectively,  is  required. 

Based  on  equations  (3-19)  and  (2-8)  it  is  found  that: 


=  -(  e  +K  )+ 


.+(_l)k  2Jjk  2U 


+( -1 )k_1 Jk~1U. 


~  2  ,  ,,  ,  sk-2  k-2T7  ,,  .  sk-1  „.k-lT7 

D  =  «>  e.+...+(-l)  e  0u  U  o+(-l)  “  U,  0 

1  1  m-2  k-5  m-1  k-2 


+  (-l)k«'kU 


k-1 


B2  =  -<ue1  +  a,2e2U2  +  *  ‘  •+(-1)k~2“k  2uk-2  +  (‘1)k~la)k  luk-l 


C2  =  -uK 


n  ,  2  r7  ,  _  .k-2  k-2TT 

D2  =  -.eQ^  eiU2+...  +  (-l)  em_2“  Uk_2 


+  (-l)k  1e  1“’k'1U  +  (-l  )k|”kU 

m-1  k-1  k 


(3-20) 


and  from  equations  (2-9)  it  is  found  that: 


P  = 


B2D1“B1D2 

C2B1 


(3-21) 


For  a  type  1  system,  eQ  in  equations  (3-20)  is  set  equal  to 
zero,  for  a  type  2  system  eo  =  e-^=0,  etc.  On  the  basis  of 
equations  (3-20)  and  (3-21)  a  cascade  compensator  can  be 
designed . 

Example  3-5 

For  the  system  of  Figure  (3-3)  it  is  desired  to  design  a 
cascade  compensator  which  places  a  pair  of  roots  at  zeta=0.5 
and  omega=l.  The  error  coefficient  Kg  should  be  50. 
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and  since  V  =  ^ ,  Z=1.529.  This  is  a  lag  network  for  which  the 
factor  0.0117  is  inherent  in  the  R-C  filter  design. 

Although  treatment  will  not  be  presented  here,  the 
algebraic  application  of  the  parameter  plane  technique  can 
be  readily  applied  to  combination  cascade  and  feedback 
compensat ion . 

B.  DOM I NANCY  OF  THE  SPECIFIED  ROOTS 

In  the  preceding  section  nothing  was  done  in  the 
calculations  to  make  the  specified  roots  a  dominant  pair. 

As  mentioned  earlier,  the  ability  to  predict  a  system's 
response  on  the  basis  of  the  location  of  a  pair  of  complex 
conjugate  roots  was  based  on  the  assumption  that  the 
magnitude  of  the  real  part  of  the  specified  roots  was  much 
less  than  that  of  the  remaining  roots  of  the  characteristic 
equation.  In  practice,  if  the  real  part  of  the  specified  or 
primary  roots  is  one  half  to  one  fifth  or  less  of  the  real 
parts  of  all  secondary  roots,  the  system  is  said  to  be 
dominant  in  the  primary  roots.  In  many  cases  the  system 
will  still  meet  the  specifications  even  if  two  pairs  of  complex 
roots  have  the  same  real  part,  provided  the  zetas  for  both 
pairs  of  roots  meet  the  specifications,  and  the  undamped 
natural  frequencies  are  such  that  the  component  time  responses 
are  not  highly  additive.  Further,  even  if  there  exists  a 
characteristic  root  whose  real  part  is  closer  to  the  origin 
than  that  of  the  primary  pair,  the  presence  of  a  closed-loop 
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zero  could  make  the  residue  of  the  close-in  root  negligible 
as  compared  to  the  primary  root.  If  possible,  however,  one 
usually  attempts  to  make  the  real  parts  of  all  secondary 
roots  large  in  magnitude. 

In  the  preceding  examples  it  should  be  pointed  out  that 
in  many  cases  there  were  actually  three  and  possibly  four 
variable  parameters.  For  instance,  the  forward  path  gain  was 
usually  a  fixed  value  in  the  computations  so  as  to  meet  the 
minimum  steady-state  accuracy  requirements.  There  is, 
however,  no  reason  why  the  gain  cannot  be  raised  above  the 
minimum  value,  thus  permitting  a  third  degree  of  freedom. 

When  cascade  and  feedback  compensation  are  used  simultaneously, 
the  forward  path  gain  and  tachometer  gain  become  the  third 
and  fourth  parameters. 

Recall  that  the  system  characteristic  equation  has  the 
m  k 

form  f(S)  =  I  akS  =0 ,  where  ak=bkQ+ck^+dk •  To  realize  the 
k — 0 

system  specifications,  one  places  a  pair  of  complex  roots 
at  S=-  *^1+ ^  ,  which  implies  that  S  +2 C ^ ^=0 . 

Since  and  are  known,  the  quadratic  can  be  divided 
out  of  the  characteristic  equation,  leaving  a  polynomial 
which  contains  the  secondary  roots  of  the  characteristic 
equation.  Since  only  two  of  the  variable  parameters  were 
used  in  fixing  the  primary  roots,  the  remaining  parameters 
will  appear  in  the  coefficients  of  the  quotient  polynomial, 
and  it  is  these  parameters  that  can  be  varied  to  achieve 
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dominance . 


Instead  of  division  to  obtain  the  quotient  polynomial, 


coefficients  of  like  powers  will  be  equated  to  achieve  a 


system  of  equations.  Let  the  quotient  polynomial  be 


given  by: 


f,(s)  =  Z  f.S*  =  0 
1  k=0  K 


(3-22) 


where  n=m-2 ,  i.e.,  equation  (3-22)  is  of  order  two  less  than 


the  characteristic  equation.  Applying  equations  (2-1), 


(3-22),  and  the  quadratic  it  follows  that: 


o  2  n  k  m  k 

(S^+2C.u.S  +  ^)(  z  f,SK)  =  £  a,  SK 

11  1  k=0  K  k=0  K 


(3-23) 


Taking  a^=l  and  equating  coefficients  of  like  power: 


a  =  f  =1 
m  n 


a  =  f  .+2?  u 
m-1  n-1  11 


a  ,,  =  f  9+22,w  f  .-Ho 
m-2  n-2  * 1  1  n-1  1 


( 3-24) 


a0  =  f  +2  f  w  f,  +f0J  . 

2*  O  ill  ^  1 


al  ■  2'-l“l1o  +  1l"l 


a0  ‘  Vl 


*  f  */  V.V,  •'  /.V.A  r  m  *  f  a  ^ .  "  V  *  *fV  VV.tf  *  ‘'a'1".  \ 
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Equations  (3-23)  and  (3-24)  can  be  solved  for  the  coefficients 
f  in  terms  of  the  coefficients  a.  The  results  will  be 


applied  to  the  following  cases: 

Case  of  k=3,  n=l 
Equation  (3-23)  becomes: 

(S2+2C1^1s+^)(f1S+f0)=S3  +  a2S2  +  a1S+a0 

Equating  coefficients  of  like  power  one  obtains: 


a3  *  fl  "  1 


a2  ’  fO+2Sl“l£l 


a  s  f  (u  +2C  co  f 

ai  X1  1  “^1  1X0 


an  in 


l0  0*1 

Solving  for  the  coefficients  f  results  in: 


f,  =  1 


f  n  n 


VP 

2?l"l 


a2-2 S^l 


(3-25) 
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Case  of  k=4 ,  n=2 


Proceeding  as  before,  the  a  coefficients  become: 


a4  f  2  1 


a3  2Cla)l  +  fl 


a2  =  V2Vlfl+fo 


(3-26) 


flWl+2Vlf0 


f  u2 
o  1 


When  solved  for  the  coefficients  f,  equations  (3-26)  yield: 


f2  =  1 


a9  “l 

a3  H  1  2C1JJ1  2CX 


^0  =  al  _  Vl  +  (J)2 

2  2  £  a  2  £  1 

i  i  n 


ao  l  ,  25iaos 

2  ,  3  2''al  oo,  ; 

11  ^  i  1 


(3-27) 


Case  of  k=5 


Similarly : 


a5  f3  1 


a4  =  2Cl("l  +  f2 


a3  ="l  +2^lVfl 


a2  =  “i  f2+24,oj  1f1+fQ 


al  =  Vl+2«l“lfO 


ao  aiifo 


a9  4|2  1  2f  a 

a^-24,-0,  =  -4  +  a„(  — -  -r) - M- 


“’1  1  ,,2  0  •  ,4 


OJ 

1  1 


2  £  ^ 

2*1  1 


al 

_  + 

ao 

“l 

2  C  a) 

3 

u.4 

2| 

1 

1 

1 

ai  2^  .  an  2  2  2 

*4 - i-£  =  a„  -  2£,o>  a  +  4C^2  -  o 
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Although  the  coef f ici ents  o f  f  have  been  derived  for  only  up  to 
the  fifth  order  case,  they  can  easily  be  obtained  for  higher 
order  cases  if  necessary. 

Example  3-6  (Third  Order  Character ist i c  Equation) 


Design  a  cascade  compensator  for  the  system  of  Figure  (3-4) 
to  obtain: 

1.  Characteristic  roots  at  zeta=0.5  and  omega=40. 

2.  K  >250. 

e— 

3.  The  specified  roots  are  to  be  dominant. 

The  characteristic  equation  of  Figure  (3-4)  is: 

S3  +  (4+P)S2  +  (4P+K?)S  +  KP  =  0 
or 

S3  +  ( 4+  °  )S2  +  (4*+K/3)S  +K«  =  0 


h""l 


where  *  =P  and  (3  =  y. 


Here  G 


=  =  30  e°=°’  6l=4>  62  =  1,  U2=1,  U3=°- 

Application  of  equations  (3-20)  yields: 

B1  =  -K+co2  =  -K+1600  B2  =  -4co  +  2Lt2  =  1440 

C  =  0  C2  =  -oo  K  =  -40K 

D1  =  4oj2  -  co3U2  =  -57600  D2  =  4u)2U2  -  o3U3  =  6400 

From  equat ions ( 3-21 )  are  obtained: 

„  _  57600 
-K+1600 

r  _  1440(-57600)-(-K+1600)( 6400 )  ,0  „ox 

-40K( -K+1600) 

For  any  value  of  K,  equations  (3-28)  will  produce  values  of 
alpha  and  beta  that  provide  characteristic  roots  at  zeta=0.5 
and  omega=40.  However,  only  certain  ranges  of  K  will  meet 
the  steady-state  error  and  dominance  requirements .  To  satisfy 
the  error  consideration  it  is  necessary  that  K>1000.  Since 
f^(s)  is  of  order  one,  equations  (3-25)  apply  and: 

f  2, 

( a, -w  ) 

2  £  1  2  11 
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From  the  characteristic  equation  it  is  seen  that 


Example  3-7  (Fourth  Order  Characteristic  Equation) 


K 

(s+oY)(s+T)(s+5')(s  +7o) 
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Figure  3-5 

Compensate  the  system  of  Figure  (3-5)  using  tachometer  plus 
acceleration  feedback  to  obtain: 

1.  Characteristic  roots  at  zeta  =  0.5  and  omega  =  2. 

2.  K  >12. 

e— 

3.  Dominance  of  specified  roots. 

The  characteristic  equation  of  Figure  (3-5)  is: 

S4+16.5S3+(73+a)S2+(82.5+B)S+25+K  =  0 


where  u=KK  and  8=KK  . 

a  t 


G  = 


;(s)  S'1  +  lG.5S3  +  73S2*82.5S  +  25 


•17 
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By  inspection  eQ=25,  e^=82.5,  e2=73,  e3=16.5,  e^=l,  U2~l, 

U3=0,  and  U4=-l.  When  equations  (3-9)  are  employed  one  finds: 

„  _  K-135  q  _  K-24 

4  y  2 

From  the  characteristic  equation  it  is  seen  that  a4=l,  a3=16.5 

a2=73+a,  a4=82.5+8,  aQ=25+K.  Since  the  quotient  polynomial 

2 

f^(s)  is  a  quadratic,  i.e.,  S  +f3S+f2=0,  equations  (3-27)  apply 

It  would  be  desirable  that,  from  a  dominancy  viewpoint, 

f  =  5.  However,  looking  at  the  dominancy  equations 

for  this  case  (equations  (3-27)),  it  is  seen  that  one  of 

the  several  expressions  for  f^  is  f^=a.j-2|j^,  which  is  a 

fixed  constant  even  though  the  remaining  expressions  for  f^ 

involve  one  or  more  variables.  Thus,  f1=14.5.  Noting  the 

a0 

most  simple  expression  for  fQ  in  equations  (3-27),  fQ=  — ?r* 

Now,  since  f1=14.5>5,  a  dominant  situation  already  exists. 
However,  the  system's  performance  can  be  further  improved 
by  choosing  appropriate  values  of  zeta  and  omega  for  the 
secondary  roots.  From  the  error  specification  it  is 
necessary  that  ^->12  or  K>300.  Now  f  ( s  )=S2  +  14 . 5Sh  ^ 

or  f  (s)=S2+14.5S+6.25+0.25K.  For  K=300,  f;[(s)=  2 

9  2 

S  +14.5S+81.25.  Therefore,  2l2^2  =  14.5,  u)2  =  81.25  implying 

^  =9.  Then  zeta=0.806.  These  appear  to  be  reasonable 

values  for  l2  and  '->9  since  the  secondary  roots  taken  aJone 

would  produce  less  overshoot  and  a  smaller  settling  time 

than  the  primary  roots.  Using  this  value  of  K  one  obtains 


a  =  41.25, 


6  =  138 


and  since  a=KK,  and  S=KK  ,  K  is  found  to  be  0.1375  while 
t  a  t 

K  =0.46. 
a 

As  an  added  bonus  of  this  method,  all  the  roots  of  the 
characteristic  equation  are  now  known  and  the  system's 
time  response  could  be  computed  if  desired. 


IV.  PROGRAM  DESCRIPTION 


The  goal  in  developing  any  computer  aided  design  program 
should  be  twofold:  (1)  provide  the  user  with  a  single, 
easily  understandable,  easily  usable  and  comprehensive 
engineering  tool,  and  (2)  dramatize  its  efficiency  above 
that  of  other  currently  available  methods.  Through  the  examples 
of  the  following  chapter  the  second  goal  will  be  demonstrated. 

It  is  first  desirable  to  reveal  the  methodology  and  internal 
structure  of  the  parameter  plane  curve  program — as  a 
consequence  it  is  hoped  that  the  first  goal  will  be  affirmed. 

A.  THE  PROGRAM 

The  parameter  plane  curve — generating  program,  or 
"program"  as  it  will  be  called  henceforth,  consists  of  a 
large  driving  routine  which  includes  all  necessary 
calculations  with  which  to  generate  the  curve  data,  and 
several  supporting  subroutines  (i.e.,  curve  plotting,  root 
solving,  data  saving,  etc.).  This  entire  package  is 
included  as  a  user-selected  option  within  another  controls 
system  computer  aided  design  package.  Among  other  options, 
the  latter  CAD  program  includes  a  root-locus  analysis-- 
as  mentioned  earlier,  the  usefulness  of  either  the  parameter 
plane  or  root  locus  technique  for  design  of  a  controls 
system  is  somewhat  limited,  but  in  combination  their 
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effectiveness  is  synergistic  (Chapter  V  will  assert  the  dual 
roles  of  the  root-locus  and  parameter  plane  methods). 
Facilities  available  within  the  program  are  many;  the 
major  options  are: 

1.  Plotting  of  constant  zeta  curves,  with  zeta  as  a 
function  of  omega. 

2.  Plotting  of  constant  omega  curves,  with  omega  as  a 
function  of  zeta. 

3.  Plotting  of  constant  sigma  (real  root)  curves. 

4.  Plotting  of  constant  zeta-omega  curves. 

5.  Tabular  output. 

6.  Rescaling  of  the  plots. 

Input  of  certain  data  is  required  to  enable  the  program. 
These  inputs  include: 

1.  Starting  value  of  . 

n 

2.  Decades  of  w  to  be  considered. 

n 

3.  Number  (and  values)  of  constant  zeta,  omega,  sigma, 
and  zeta-omega  curves. 

4.  Coefficients  associated  with  the  constant,  alpha,  and 
beta  terms  of  the  characteristic  equation. 

Each  of  the  basic  program  option  areas  will  be  described  in 

appropriate  detail,  as  well  as  their  interaction  with  the 

input  data. 

1 .  Constant  Zeta  Contours 

In  practice  design  specifications  for  control 
systems  are  given  in  terms  of  percent  overshoot,  settling 


time,  error  constraints,  etc.  A  value  of  zeta  can  be 
associated  with  the  first  of  these  spec i f icat ions--t hat  is, 


m 


a  given  percent  overshoot  requirement  can  be  related  to  a 
specific  value  of  zeta.  Given  a  specific  zeta  value,  the 
program  calculates  the  alpha  and  beta  coefficients  of  the 
characteristic  equation  by  holding  the  zeta  value  constant 
and  varying  the  value  of  omega.  The  limits  within  which 
omega  is  varied  are  defined  by  the  user's  choice  of  the 
initial  value,  and  the  number  of  decades  of  omega  to 
be  considered. 

From  the  nature  of  the  mapping  process,  it  is  clear 
that  when  the  contour  of  the  coefficient  plane  passes  through 
a  designated  point  (M-point),  the  original  mapping  contour 
on  the  S-plane  passes  through  a  point  which  is  a  root  of  the 
characteristic  equation.  The  zeta  value  chosen  for  the 
contour  is  then  the  zeta  for  the  root.  The  value  of  omega 
associated  with  the  M-point  is  the  radial  distance  from  the 
origin  of  the  S-plane  to  the  root.  Thus,  a  complex  root 
is  determined  when  the  M-point  lies  on  a  constant  zeta  curve 
of  the  parameter  plane.  The  value  of  this  root  and  its 
complex  conjugate  is: 

S  =  -£W+ja>/  l-£2 

If  the  characteristic  equation  is  such  that  several  complex 
roots  exist,  then  the  parameter  plane  curves  required  to 
realize  these  roots  must  all  pass  through  the  M-point.  If 

i 

the  complex  roots  have  the  same  zeta  but  different  omega 


values,  then  the  constant  zeta  curve  must  pass  through  the 


M-point  more  than  once.  In  fact,  once  any  point  on  the 
zeta  contour  is  defined,  omega,  alpha,  and  beta  are  also 
defined,  and  all  roots  of  the  characteristic  equation  are 
thus  fixed. 

2 .  Constant  Omega  Contours 

Within  the  program,  for  a  given  value  of  omega,  zeta 
is  varied  between  zero  and  one  inclusively  while  omega  is 
held  constant. 

As  with  the  constant  zeta  curves,  any  point  on  a 
constant  omega  contour  is  the  omega  for  a  complex  root  of  the 
characteristic  equation.  By  selecting  an  operating  point, 
zeta  is  also  defined  whereby  a  pair  of  complex  conjugate  roots 
is  established.  Again,  once  any  point  is  chosen  on  either 
a  constant  zeta  or  a  constant  omega  curve,  all  roots  of  the 
characteristic  equation  are  established. 

3 .  Constant  Sigma  Contours 

When  real  roots  are  to  be  evaluated,  it  is 
convenient  to  return  to  the  characteristic  equation: 

m  k 

2  a,  SK  =  0  (2-1) 

k=0  K 


where,  again,  a^  represents  a  linear  combination  of  constant, 
alpha,  and  beta  terms.  S=-<r  (a  real  number)  is  then  an 
equation  of  a  straight  line  on  the  parameter  plane.  If  any 
line  of  constant  sigma  value  passes  through  the  M-point, 


then  the  alpha  and  beta  coordinates  of  the  M-point  satisfy 


the  characteristic  equation  for  a  real  root  located  at  -<r. 

For  program  considerations,  one  enters  a  positive 
value  for  sigma  (corresponding  to  a  real  root  at  -<?)  and 
the  constant  sigma  contour  (straight  line)  is  developed. 

The  coordinates  of  any  point  on  this  curve  produce  a  real 
root  at  -<r.  That  this  is  a  useful  tool,  consider  that 
system  specifications  can  be  achieved  by  placing  a  pair  of 
complex  conjugate  roots  of  the  characteristic  equation  at 
a  specific  location.  To  ensure  dominance  of  this  root 
pair,  the  real  part  of  the  complex  roots  so  placed  should 
be  smaller  in  magnitude  than  that  of  the  remaining  system 
roots.  Roots  placed  at  a  specified  sigma  value  can  thus 
ensure  at  least  one  real  root  whose  magnitude  is  greater 
than  the  real  part  of  the  intended  complex  conjugate  pair. 

4 .  Constant  Zeta-Omega  Curves 

For  a  fixed  value  of  zeta  and  omega  a  pair  of  complex 
conjugate  roots  is  defined  in  terms  of  the  expression: 

2 

S  =  +  V  1-c, 

The  real  part  of  these  roots  is,  thus,  defined  by  the 

zeta-omega  product.  Note  that  settling  time  is  defined  as 
4 

T  =  —  .  If  the  product  is  known,  so,  too,  is  the 

s  |gj 

duration  of  the  transient  response.  Thus,  by  specifying  a 


constant  value  for  the  zeta-omega  product,  any  point  on  the 
contour  generated  by  this  value  will  produce  a  given 
settling  time. 


For  any  of  the  parameter  plane  contours,  it  is  desirable 
to  ascertain  the  values  of  the  characteristic  roots  for  every 
few  values  of  alpha  and  beta.  This  feature  has  been 
incorporated  within  the  tabular  output  facility  as  described 
below. 

5 .  Tabular  Output 

For  each  of  the  zeta,  omega,  and  zeta-omega  parameter 
plane  curves,  an  arbitrary  though  reasonable  300  points  are 
calculated  with  which  to  plot  the  contours.  For  the  constant 
sigma  curves,  only  two  points  are  needed  to  define  the 
required  straight  lines  (in  practice,  4  points  are  generated 
to  ensure  that  the  sigma  contours  can  be  plotted  within  the 
user-defined  axes  limits).  Because  of  the  bulk  of  data 
points  so  generated,  tabular  output  is  offered  as  an  option 
(as  is  graphical  output),  and  all  points  so  generated  are 
listed  for  the  user.  In  addition,  it  is  worthwhile  to 
calculate  system  roots  for  given  values  of  alpha  and  beta. 
However,  computation  of  roots  for  each  alpha  and  beta  pair 
would  cost  unnecessary  computer  time  and  will  likely  tax  the 
user's  patience  with  the  bulk  of  output  so  generated.  Thus 


the  characteristic  roots  are  generated  for  every  tenth 
pair  of  alpha  and  beta  values. 

6 .  Plot  Rescaling 

Regardless  of  the  plot  scale  selected  by  the  user, 
the  program  generates  the  full  300  data  points  for  each 
curve  requested  (4  points  for  constant  sigma  lines).  When 
the  graphical  output  option  is  requested,  the  first  family 
of  curves  is  automatically  scaled  to  encompass  each  and 
every  data  point.  The  disadvantage  of  this  technique  is 
that,  because  most  activity  for  the  vast  majority  of  systems 
occurs  near  the  physical  origin  (i.e.,  alpha=beta=0) ,  the 
curves  may  at  first  appear  within  only  a  very  small  sector 
of  the  entire  plot  area,  and  often  they  are  indistinguishable 
from  one  another.  The  advantages  of  automatic  scaling  for 
the  first  set  of  parameter  plane  curves  far  outweigh  this 
disadvantage.  First,  by  plotting  all  available  data  points, 
the  possible  limits  for  alpha  and  beta  are  exposed — this  is 
important  if  very  large  values  of  alpha  and/or  beta  are 
required  to  meet  the  design  specifications.  Second,  for 
some  systems  the  area  of  activity  may  not  occur  near  the 
origin,  and  automatic  scaling  spares  the  user  the  task  of 
selecting  a  sector  and  possibly  missing  a  sector  of  interest. 


Although  a  seemingly  arbitrary  choice,  the  generation 
of  characteristic  roots  for  every  tenth  alpha,  beta  pair 
produces  a  very  tidy  output  on  the  commonly-used  IBM-3278 
computer  terminal. 


Once  the  user  is  able  to  view  the  panoramic  alpha, 
beta  parameter  curves,  it  becomes  obvious  which  sector(s) 
are  of  interest.  The  user  then  has  the  option  to  rescale 
the  set  of  curves  by  selecting  upper  and  lower  limits  for 
the  alpha  and  beta  axes.  He  may  continue  to  rescale  the 
family  of  parameter  curves  as  often  as  is  desirable,  and  at 
any  time  the  autoscaling  option  may  be  recalled. 

The  curves  generated  by  the  above  program  are 
sufficient  to  explore  most  control  system  engineering 
problems.  The  use  of  these  parameter  plane  contours,  and 
their  interaction  with  one  another,  will  be  evidenced  in 
the  next  chapter.  The  source  code  listing  of  the  program 
is  included  as  Appendix  B. 


B.  INSTRUCTION  TO  THE  USER 

The  parameter  plane  program  is  highly  interactive--the 
user  is  prompted  for  each  required  input.  A  brief 
description  of  all  but  the  most  trivial  input  items  follows. 


-  Starting  value  of  &  :  For  most  control  systems  the 
initial  value  of  <*>  is  chosen  to  be  zero.  Because  u 
is  used  in  the  denominator  of  certain  of  the  paramet§r 
plane  equations,  ^  must  be  greater  than  zero.  However, 
the  user  may  choos§  w  arbitrarily  close  to  zero  if 
desired. 


-  Number  of  decades:  For  the  majority  of  control  system 
problems  a  suitable  number  of  decades  to  be  considered 
might  be  two  or  three.  For  higher  order  systems,  it 
would  be  advisable  to  start  with  a  slightly  larger 
number  of  decades,  especially  if  the  initial  ^  value 
is  small.  For  subsequent  families  of  curves,  Phe 
number  of  decades  can  easily  be  changed. 


-  Constant,  alpha,  and  beta  coefficient  values: 
Characteristic  coefficients  are  requested  from  the 
highest  to  lowest  order  term.  By  way  of  an  example, 
a  third  order  characterist ic  equation  might  be: 

S3  +  (3a  +  3  +  10)S2  +  as  +  (3+5)  =  0 

Here,  the  constant  coefficients  would  be  entered  in  the 
following  sequence:  1,10,0,5  while  alpha  and  beta 
coefficients  would  be  entered  as  0,3, 1,0  and  0,1, 0,1 
respectively . 

-  Zeta  values:  By  convention,  values  are  restricted  to 
between  zero  and  one,  inclusively. 

-  Sigma  values:  Positive  values  of  sigma  correspond  to 
negative  real  roots.  Since  few,  if  any,  practical 
engineering  applications  exist  for  designing  a  positive 
real  root  into  a  system,  negative  values  for  sigma  are 
disallowed . 

-  Omega  values:  Values  for  constant  omega  curves  are 

restricted  in  the  lower  limit  to  the  starting  +  value, 
and  in  the  upper  bound  by  ojnxl0°eca°es  .  n 

-  Zeta-omega  values:  As  with  the  constant  sigma  curves, 
values  for  constant  zeta-omega  contours  must  be  greater 
than  or  equal  to  zero. 

The  user  then  has  the  following  options: 

1.  Review  entries. 

2.  Change  any  entry. 

3.  Tabular  output. 

4.  Graphical  output. 

Remember  that  tabular  output  includes  300  data  points 
for  all  but  the  sigma  contours.  Characteristic  roots  are 
displayed  for  every  tenth  alpha,  beta  pair.  Because  of  the 
bulk  of  output  for  this  option,  use  it  only  when  necessary. 
If  a  printed  copy  of  the  tabular  output  is  desired,  type 


"record  on"  before  invoking  the  program.  Upon  exiting  the 
program,  type  "record  off",  after  which  the  user  may  save 
the  preceding  terminal  session  in  a  listing  file  designated 
by  a  name  of  his  choosing.  Simply  print  the  listing  file, 
which  will  include  all  output  which  has  transpired  on  the 
terminal  between  the  two  calls  to  "record". 

When  graphical  output  is  requested,  all  curves  are 
superimposed  on  the  same  plot.  The  first  set  of  curves  is 
produced  with  an  autoscaling  feature,  which  plots  all 
points  calculated  (the  range  of  points  depends  on  your 
choice  of  initial  ,  number  of  decades,  zeta  values,  omega 
values,  etc.).  For  most  characteristic  equations  only  the 
first  quadrant  (i.e.,  positive  alpha  and  beta  values)  will 
be  of  interest,  since  negative  values  usually  (but  not 
necessarily)  imply  negative  characteristic  coefficients  and, 
thus,  positive  roots  leading  to  system  instability.  The 
nature  of  the  first  (autoscaled)  set  of  curves  will  reveal 
the  actual  areas  of  interest  for  subsequent  plots  for  the 
same  system. 

Finally,  after  each  family  of  curves  is  plotted,  the 
user  has  six  additional  options: 

1.  New  problem. 

2.  Same  problem. 

3.  Root  finder. 


4.  Save  problem. 


'  V.  v.  -".  ''.  ^  v.-wnv  f.v\  . 


5.  Create  "DISSPI.A  metafile". 

6.  Return  to  main  menu. 

Items  1,  4,  and  6  are  self-explanatory.  The  remaining 
options  deserve  some  additional  explanation. 

-  Option  2:  This  option  can  be  used  to  re-enter  the 
problem  at  a  point  prior  to  actual  plotting.  Then, 
specific  input  values  can  be  added  or  modified, 
tabular  and/or  graphical  output  can  be  requested,  and 
entries  can  be  reviewed.  It  is  within  this  option  that 
the  graph  coordinate  axes  can  be  rescaled.  If  user- 
defined  scaling  is  desired,  the  minimum  and  maximum 
values  for  the  axes  are  requested. 

-  Option  3:  Although  within  the  tabular  output  feature  a 
set  of  characteristic  roots  is  produced  for  every  tenth 
pair  of  alpha  and  beta  values,  the  bulk  of  output  using 
that  option  may  prove  excessive  for  some  applications. 
Here,  the  user  has  the  option  of  choosing  specific 
values  of  alpha  and  beta  (e.g.  extracted  from  the 
family  of  parameter  plane  curves)  and  obtaining  the 
system  roots. 

-  Option  5:  The  program  provides  a  choice  from  among 
four  graphic  output  devices.  Usually  the  user  will 
nominate  the  TEK618  graphics  terminal  due  to  its 
relatively  high  quality  plot  resolution.  Once  the 
user  has  produced  a  parameter  plane  plot  to  his  liking, 
he  may  wish  a  final  plot  of  very  high  resolution.  By 
selecting  this  option,  the  plot  is  stored  as  a  DISSPLA 
metafile,  and  the  program  is  terminated  ( terminat ion 

of  the  program  is  necessary  at  this  point  due  to  an 
anomaly  of  the  DISSPLA  graphics  package).  Simply  type 
"DISSPOP"  and  follow  the  simple  instructions,  choosing 
the  default  options  as  they  are  presented.  Within  the 
"DISSPOP"  routine,  any  of  several  output  devices  can  be 
called,  including  the  high  resolution  Versatec  plotter 
and  the  3800  laser  printer. 


GO 


V.  PARAMETER  PLANE  CURVES-GRAPH ICAL  METHOD 


A.  GRAPHICAL  SOLUTION 

The  algebraic  solutions  discussed  in  Chapter  III  have  the 
disadvantage  that  a  fixed  value  of  zeta  and  omega  must  first 
be  chosen  to  compute  the  alpha  and  beta  terms.  In  some 
instances  it  is  possible  to  modify  the  remainder  polynomial 
so  as  to  ensure  that  the  specified  roots  are  dominant. 

However,  it  is  not  always  possible  to  guarantee  that  roots 
placed  at  a  specified  location  can  be  made  dominant.  Thus, 
an  exhaustive  trial-and-error  procedure  may  be  required  to 
achieve  the  best  values  for  the  various  parameters.  Trial- 
and-error  may  also  be  required  in  the  design  of  cascade 
compensators  where  a  specific  root  location  may  require 
parameter  values  that  are  not  physically  realizable.  In 
these  cases,  the  calculation  must  be  repeated  in  terms  of 
slightly  modified  specifications;  possibly  a  different  means 
of  compensation  must  be  used. 

To  avoid  this  trial-and-error  analytical  approach  ,  one 
can  employ  the  graphical  solution.  Once  a  family  of  curves 
is  generated  by  the  program  one  can,  by  choosing  an  M-point 
in  the  parameter  plane,  obtain  from  the  curves  the  n  roots 
of  the  nth  order  characteristic  equation.  The  trial-and-error 
procedure  can  then  be  done  visually  to  reveal  an  operating 
point  which  best  meets  the  given  specifications. 


Example  5.1  (An  Attitude  Control  System  for  Large  Launch 


Figure  5-1 

Figure  5-1  shows  the  mechanization  for  a  control  system  fur  a 
large  launch  vehicle;  Figure  5-2  shows  the  equivalent  block 
diagram 


Figure  5-2 


where  the  equivalent  G(s)  is: 


( a„S2  +  K  s  +  ) ( S2-C ) 


S  (— 2  +  e  +  S  +  a  C) 

u  a»  l 

e  e 


From  G(s)  one  obtains  the  system's  characteristic  equation: 


6  2£ 

+  ^  S5  +  (l+an)S4  +  (a,C+K, )S3  +  (K2-aQC)S2 


U)  p 

e 


i  i 


-  CKjS  -  CK2  =  0 


where  a^ ,  a^ ,  ,  K 2  =  control  system  gains 

C  e  =  damping  ratio  of  control  servo 


“  e  =  natural  frequency  of  control  servo 


Gains  and  K 2  are  chosen  as  the  system  parameters  a  and  5, 
respectively,  to  be  portrayed  in  the  parameter  plane.  One 
must  then  find  suitable  values  for  these  parameters  to  yield 
a  desired  stability  margin  and  a  satisfactory  transient 
response.  For  a  typical  choice  of  system  parameters  (such 
as  those  describing  Saturn  V),  it  is  assumed  that: 

S  =  0.717 

e 

^  =  4.71  Hz  =  29.594  radians 


Then  the  system's  characteristic  equation  becomes: 

0.0011S6  +  0.0485S5  +  0.5S4  +  (0.7+a)S3  +  (0.35  +  S)S2 
-  0.7<*S  -  0.78  =  0 

Various  values  of  a^  and  a^  were  used  to  deduce  their  effect 
on  the  stability  regions,  which  is  indicated  in  Figure  5-3. 
The  numbers  of  stable  and  unstable  roots,  respectively,  are 
portrayed  in  parentheses  for  each  region  of  the  parameter 
plane . 

The  analysis  procedes  as  follows:  the  5=0  contour  repre¬ 
sents  boundaries  of  stability  (or  relative  stability  when 
5>0)  associated  with  pairs  of  complex  conjugate  roots.  The 
S=0  (sigma=0)  curve  represents  real  root  stability  boundaries 
The  region  of  stability  is  thus  that  area  bounded  by  these 
two  contours.  See  Figures  5-4a  and  5-4b  for  a  magnified 
view  of  this  area.  Any  negative  alpha-beta  pair  from  within 
this  region  will  exhibit  six  stable  (i.e.,  all  within  left- 
half  S-plane)  roots  and  system  stability  will  be  assured. 

Note  that  from  the  form  of  the  characteristic  equation, 
both  alpha  and  beta  must  be  negative  to  obtain  a  stable 
system.  To  illustrate,  let  us  select  an  arbitrary  operating 
point  constrained  to  lie  within  the  lower  loop  of  Figure  5-4b 


ATTITUDE  CONTROL  SYSTEM 


o 


*r  v  /.  a  ^ 


and  also  satisfying  the  requirement  that  both  alpha  and  beta 
be  negative.  Superimposed  on  Figure  5-4b  is  the  constant 
5=0.5  contour,  upon  which  our  M-point  might  be  chosen.  If 
we  select,  say,  a=-0.15  and  B =— 0 .02  ( corresponding  to  £=0.5 
and  w=0.5),  the  system  roots  are: 

-0.26  +  j  0 . 45 

-dominant  roots 

-0.26  -  j0.45 
-0.32  +  j  0 . 12 
-0.32  -  j 0 . 12 
-14.6  +  jO 
-26.7  +  jO 

Coincidentally,  the  roots  associated  with  our  choice  of  zeta 
and  omega  are  seen  to  be  dominant.  The  actual  choice  of  an 
operating  point  may  depend  on  other  criteria  not  discussed 
here . 

Since  and  K 2  (a  and  3  ,  respectively)  are  functions 
of  u  ,  5  ,  aQ ,  a^ ,  C ,  0)  ,  and  £  e ,  they  may  be  determined  for 

various  instances  of  flight  by  plotting  several  constant 
zeta  and  constant  omega  curves.  Actually,  and  K ^  vary  so 
little  within  the  range  of  values  used  for  C,  for  specified 
values  of  5  and  jJ,  that  it  becomes  possible  to  choose  constant 
values  for  and  K 2« 

This  has  been  a  relatively  simplistic  treatment  of  a 
complicated  control  system  problem,  but  it  demonstrates  the 
power  of  the  parameter  plane  graphical  te  unique.  Knowing 


Akv\ 


68 


nothing  more  than  the  system's  characteristic  equation,  the 
static  and  dynamic  stability  boundaries  may  be  obtained  to 
define  the  area  of  overall  stability.  Of  course,  other 
system  constraints  may  exist  to  further  limit  this  area. 

As  a  note,  the  root-finding  option  available  within  the 
program  was  used  to  confirm  the  numbers  of  stable  and 
unstable  roots  lying  within  each  region  of  Figures  5-3  and  5-4. 

Example  5-2  (Alternator  Voltage  Regulator) 

In  designing  a  voltage  regulator  of  the  type  shown  in  Figure 
5-5,  we  must  find  values  for  K^,  ,  and  K9  that  provide 

stability  and  good  transient  performance. 


EXCITER  ALTERNATOR 


Figure  5-5 

The  characteristic  equation  of  the  system,  including  alternator, 
tie-line,  etc.  is  (see  reference  8  for  details): 

0.0095S5  +  0.1325S4  +  1.72S3  +  (  K.,  +  7 . 55  )S2  +  ( IC  +9.1)S 

**  J. 


We  will  consider  Kq=0  as  fixed  and  a  and  3  as  variable  para¬ 
meters  corresponding  to  K2+7.55  and  K^+9.1,  respectively. 
Curves  of  a  versus  3  are  plotted  in  Figure  5-6  with  cjn  as 
the  variable  parameter  (for  this  system,  Kn  and  K2  are 
known  to  be  functions  of  wn  alone).  Since  the  program  plots 
constant  zeta  curves  as  a  function  of  varying  omega,  these 
curves  were  selected  as  logical  candidates  to  study  the 
problem. 

The  beta  axis  corresponds  to  the  zero  value  of  damping 
constant  (sigma)  since  below  the  zeta=1.0  curve,  the  real 
roots  (Thaler  and  Brown  1960)  are  the  negative  slopes  of  the 
tangents  drawn  from  the  point  in  question  to  the  5= 1  curve. 
The  machine  is  then  stable  for  any  (a,  3)  pair  between  the 
5=0  and  S=0  curves.  The  greatest  stability,  of  the  machine 
is  then  possible  when  both  zeta  and  sigma  are  largest. 
Further,  the  best  stability  can  be  expected  in  the  region 
bounded  by  the  5=0.3  and  5=1.0  contours  (Kabriel  1967). 
Similar  stability  limits  of  a  and  3  can  be  investigated  by 
taking  as  10,  20,  30  and  so  on. 

Consider  now  K^=0  as  fixed.  Then  ot  and  3  will  represent 
the  variable  parameters  K2+7.55  and  Kq-2.5,  respectively. 

The  characteristic  equation  then  becomes: 

0.0095S5  +  0.1325S4  +  1.72S3  +  aS2  +  9. IS  +3=0 


'J>  •>  -Jf  ■>  ">  ">  ->  ■>  '.-  “ ’>  -  »  '>  'J/  *.V  *T  ■>  ■>  »JI  >U(  ->  ->  -_x  ->  ->  /  w  v  v.v.v  r  -r ,  «r.  > 


Again,  Kq  and  are  known  to  be  functions  of  u)n  alone.  Here, 


however,  can  be  solved  explicitly  [Ref.  8]  and  one  obtains 


three  straight-line  equations: 


3  =  0 


5.421a  -8  =  3.894 
175.63a  -3  =  4085.0 


or  in  terms  of  system  parameters: 


Kq  =  2.5 


5 . 421( K0+7 . 55 ) -K  =  3.894  -  2.5 
2  o 


175.63(K2+7.55)-KQ  =  4085.0  -  2.5 


These  are  plotted  in  Figure  5-7.  The  triangular  region 
bounded  by  these  three  lines  represents  a  stable  region.  The 
values  of  (a ,  8 )  within  the  triangle  and  hence  corresponding 
values  of  KQ  and  K2  can  be  predicted  for  stable  operation. 

The  procedure  can  be  repeated  for  further  investigation  by 


taking  K^=15,  30,  45,  etc. 


The  analytical  parameter  plane  technique  can  be  used  to 
determine  the  stability  limits  of  Kq  and  by  choosing  several 


constant  values  of  K2 .  But  since  K2  has  to  be  selected 


arbitrarily  for  this  purpose,  this  method  becomes  cumbersome 
and  time-consuming.  The  method  presented  here,  on  the  other 
hand,  can  be  used  to  determine  the  stable  range  of  either 
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ALTERNATOR  VOLTAGE  REGULATOR 


Figure  5-6 

,  Parameter  Plane  Curves  for 
0.0095S  +  0.1325S4  +  1.72S3  +  aS2  +  BS 
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Kq  and  K0  or  and  K2  fixing  the  third  parameter.  Over  and 


above  the  ability  to  predict  -stable  operation,  the  method 
povides  a  direct  measure  of  the  damping  at  and  around  a 
chosen  operating  point. 

Example  5-3  (Lead  Compensator) 

Consider  the  system: 


K(S+Z) 

t 

S+P 

** 

S(S+1)(S-M) 

The  character ist ic  equation  is: 

S4  +  ( 5+P )S3  +  (4+5P)S2  +  ( 4P+K ) S  +  KZ  =  0 

We  choose  to  cancel  the  pole  at  S=-l  with  the  zero;  thus 
Z= 1 . 0  and  the  characteristic  equation  becomes: 

S4  +  (5  +  P)S3  +  (4  +  5P)S2  +  ( 4P  +  K  )S  +  K  =  0 

Let  P  =a  and  K  =  B.  Then  the  parameter  plane  curves  are  as 
shewn  in  Figure  5.8.  For  the  coefficients  of  the  character¬ 
istic  equation  to  remain  positive  (and  thus  ensure  stability), 
it  is  convenient  to  consider  only  positive  values  for  alpha 
and  beta.  Dy  inspection  we  choose  ^=0.5  and  u  =2.0  as  a 
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Figure  5-9 

Parameter  Plane  Curves  for 
S4  +  (5+ct)S3  +  (4+5a)s2  +  (4^+3  )s  + 3  =  o 
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good  operating  point,  for  which  a  =  4.0  and  3  =24.0.  The 
corresponding  vicinity  of  Figure  5-8  is  re-scaled  in  Figure 
5-9.  Again,  using  the  root-finding  option,  the  roots 
associated  with  this  (C  ,  ^n)  pair  are  shown  to  be  dominant. 

Example  5-4  (Lag  Compensator) 

If  we  are  especially  concerned  with  steady-state  accuracy 
for  a  ramp  input,  it  may  be  advisable  to  design  a  lag 
compensator.  The  parameter  plane  permits  us  to  consider 
steady-state  accuracy  while  designing  the  transient  response. 
If  our  system  is  the  same  as  that  considered  for  the  lead 
compensator,  the  characteristic  equation  remains: 

S4  +  (5+P)S3  +  (4+5P)S2  +  ( 4P+K )S  +KZ  =  0 


But  now  the  error  coefficient  is  Kg  =  ^  .  Having  three 
unknown  parameters,  K,  Z,  and  P,  two  must  be  selected  (or 
some  combination  of  two)  to  be  a  and  8  while  a  numerical  value 
is  assigned  to  the  third.  Once  this  choice  has  been  made 
thr'  parameter  plane  curves  can  be  calculated  and  plotted,  and 
the  loci  of  constant  can  be  superimposed.  Let  Z=0.1, 

P=  a,  and  K=  8.  The  characteristic  equation  becomes: 

S4  +  (5+°0S3  +  ( 4  +  5  &)S2  +  ( 4  cm-  3)S  +  O.lB  =  0 


m 


Parameter  plane  curves  are  shown  in  Figures  5-10  and  5-11. 

0  1  3 

Lines  of  Kg  =  =  0.1,  0.2,  0.3,  etc.  may  be  superimposed. 

If  we  select  K&  =  0.15  and  (similarly  to  the  lead  compensator 
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example)  5=0.5  and  o>  =2.0,  then  a  =2 . 14  while  3=12.89.  This 

n 

produces  dominant  roots  at  -1.0  +  jl.7. 

Let  us  reconsider  the  systems  for  which  a  lead  compensator 

and  a  lag  compensator  were  designed.  The  characteristic 

equation  as  well  as  the  error  coefficient  each  contain  three 

unknowns  (if  we  assume  a  numerical  value  for  K  ).  We  can 

e 

imbed  the  error  coefficient  in  the  characteristic  equation 
by  direct  substitution,  thereby  eliminating  one  of  the 
unknowns.  Let  us  eliminate  the  gain  parameter  K  -  note  that 


4PK 


Returning  to  the  characteristic  equation: 

48  p  4PK 

S  +  (5+P)S  +  (4+5P)S  +  (4P+-b-£)S  +  4PK  =  0 

&  e 

p 

Letting  P  =  a  and  ^  =  3  ,  the  characteristic  equation  becomes: 

S4  +  ( 5-kx  )S3  +  (4+5a)S2  +  (4a+4KeB)S  +  4K  a  =  0 

If  a  constant  value  is  now  chosen  for  Kg ,  say  Kg  =  2.0,  the 
characteristic  equation  becomes: 

S4  +  C  5+ct  ) s 3  +  (4+31  )S2  +  (4a  +8  3)S  +  8a  =0 

for  which  the  parameter  plane  contours  are  first  shown  in 
Figure  5-12  and  are  further  magnified  in  Figure  5-13.  Now 
when  any  operating  point  is  chosen  on  the  parameter  plane 
curve(s),  the  selected  (a  ,3  )  pair  generates  Ke=2.0.  Of 
course,  this  procedure  can  be  repeated  for  any  choice  of  K  . 


Figure  5-13 

Parameter  Plane  Curves  for 
( 5  +  a)S3  +  (4  +  50OS2  +  (4a+88)S  +  8a 


For  this  system,  the  following  specifications  are  to  be  met: 


1.  Set  K  to  the  stability  limit. 

2.  Place  a  dominant  root  pair  within  the  following  region: 

0 . 4<C  ^_0 . 7  ,  and  2<w <_6 . 

3.  Both  tachometer  and  acceleration  feedback  may  be  used, 
but  if  possible  choose  only  one. 

From  the  figure  the  uncompensated  system's  open  loop  transfer 

function  is: 

GH  =  -1  =  - ^ - 

S(S+10)(S  +5S+100) 

From  which,  when  expanded,  the  characteristic  equation 
becomes : 

S4  +  15S3  +  150S2  +  1000S  +  K  =  0 

To  determine  the  value  of  K  at  the  stability  limit  the  Pouth 
array  is  employed: 


1 

150 

K 

15 

1000 

0 

1250 

15K 

0 

1.25x106-225K 

0 

0 

Here,  the  stability  limit  is  seen  to  be  K=5555.5.  If  both 
tachometer  and  acceleration  feedback  are  used  the  compensated 
system's  characteristic  equation  becomes: 

S4  +  15S3  +  (150+5555.5a)S2  +  ( 1000+5555 . 53 ) S  +  5555.5=0 

where  a=K  ,  g=K  ,  and  K  has  been  set  to  the  stability  limit. 
Parameter  plane  curves  for  this  system  are  shown  in  Figure 
5-14. 

From  these  curves,  the  following  analysis  can  be  made. 

The  origin  of  the  parameter  plane  corresponds  to  the  roots 

of  the  uncompensated  system.  Since  the  £=0  curve  passes 

through  the  origin,  two  roots  are  located  on  the  ja>  axis  of 

the  S-plane  as  was  to  be  expected  from  the  Routh  array.  The 

remaining  two  roots  are  also  complex  and  correspond  to  £=0.8 

and  co=5.0.  It  is  important  to  note  that  when  an  operating 

point  involves  two  different  pairs  of  complex  roots,  then  the 

curves  for  two  different  values  of  omega  and  two  different 

values  of  zeta  must  pass  through  the  point.  To  determine 

which  value  of  omega  corresponds  to  which  value  of  zeta,  it 

becomes  necessary  to  refer  to  the  program's  tabular  data 

output,  which  is  not  included  here  due  to  lack  of  space. 

With  K  =0,  the  effect  of  tachometer  feedback  alone 
a 

corresponds  to  movement  of  the  M  point  along  the  8  axis. 

In  Figure  5-14  the  unstable  region  is  determined  by  an 
inspection  of  the  way  the  constant  zeta  curves  tend 
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as  zeta  increases.  Since  the  6  axis  is  always  in  the 

unstable  region,  it  is  concluded  that  tachometer  feedback 

alone  cannot  stabilize  the  system. 

The  effect  of  acceleration  feedback  alone  can  be  observed 

by  traveling  along  the  a-axis  of  Figure  5-14.  If  K  is 

£1 

varied  between  0.01  and  0.06,  the  system  will  exhibit  two 
pairs  of  complex  roots  with  the  following  ranges  of  values 
for  zeta: 

0. 3<£  <0.5  and  0.25<£<0.32 

If  both  tachometer  and  acceleration  feedback  are  used, 

it  is  seen  that  tachometer  feedback  will  in  general  cause 

the  zeta  of  one  pair  of  roots  to  increase  while  the  zeta 

of  the  other  pair  decreases.  Acceleration  feedback  alone 

would  appear  to  be  the  better  choice. 

From  the  set  of  curves  it  is  determined  that  with  K  =0 

t 

and  Ka=0.012,  four  complex  roots  are  located  with  associated 
values  at  £=0.45,  <jj=4.0,  and  £=0.32,  w=13.0.  Since 
( 0. 45  )  ( 4  )  =  L8< < ( 0 . 32  )  ( 13 )  =  4 . 15  ,  it  is  apparent  that  the  roots 
at  £=0.45  and  w=4.o  are  dominant.  The  specifications  have 
been  met  and  the  problem  is  solved. 

B.  MISCELLANEOUS  ASPECTS  OF  THE  PARAMETER  PLANE 

It  can  be  demonstrated  that  constant  zeta  parameter 
plane  curves  of  order  two  through  five  originate  at  a  point 
where  alpha  =  M  and  beta  =  where  M,  N,  and  K  are 

K.  /V 


determined  by  the  zero  and  first  power  coefficients  only. 
Intuition  can  be  used  to  conclude  that  constant  zeta  curves 
of  any  order  originate  at  this  common  point  which  is 
determined  only  by  the  zero  and  first  power  coefficients. 

An  exception  is  when  K=0.  In  this  case  the  origin  point 
depends  on  higher  order  coefficients  and  its  location  will 
be  obvious  given  a  specific  problem.  If  K  is  not  zero, 
the  origin  is  independent  of  the  order  of  the  characteristic 
equation. 

Inspection  of  the  expressions  for  alpha  and  beta  indicates 
that  the  shape  of  the  constant  zeta  curves  as  omega  becomes 
larger  is  primarily  determined  by  the  coefficients  of  higher 
power,  and  in  general  the  curves  become  more  complex  and 
less  well  behaved  as  the  order  of  the  characteristic 
equation  increases.  For  a  given  characteristic  equation, 
an  increase  in  complexity  can  be  observed  as  alpha  and  beta 
appear  in  more  coefficients. 

All  constant  zeta  curves  tend  to  plus  or  minus  infinity. 
The  relative  magnitudes  of  the  coefficients  determine  whether 
the  limit  is  plus  or  minus  infinity.  It  is  therefore 
necessary  to  choose  a  frequency  range  of  interest  before 
plotting  the  curves,  thus  limiting  the  analysis  to  one 
"window"  of  the  infinite  plane. 

Since  no  stability  criteria,  either  relative  or 
absolute,  has  been  established  for  the  parameter  plane,  it 


is  necessary  to  base  the  stability  analysis  on  observing 
which  way  the  curves  tend  as  omega  and  zeta  are  varied. 

For  this  reason  it  is  worthwhile  to  plot  curves  for  as  many 
values  of  zeta,  omega,  sigma,  and  if  desired,  zeta-omega, 
as  are  necessary  to  ascertain  the  pattern. 


VI.  CONCLUSIONS 


Parameter  plane  techniques  have  been  applied  to  the 
compensation  of  linear  control  systems.  General  equations 
have  been  derived  for  the  cases  of  feedback,  cascade,  and 
combination  feedback-cascade  compensation,  to  enable  one  to 
place  a  pair  of  complex  conjugate  roots  at  a  specific 
location  in  the  S-plane,  while  simultaneously  satisfying 
the  steady-state  accuracy  requirements.  A  dominancy 
technique  has  been  introduced  whereby  once  a  pair  of  complex 
roots  is  fixed,  the  remaining  roots  of  the  characteristic 
equation  can  be  manipulated  to  ensure  that  the  specified 
roots  are  dominant. 

The  impetus  for  development  of  a  parameter  plane  program 
was  to  provide  the  user  with  a  quick,  simple  means  of 
obtaining  the  information  available  in  the  analytical 
solution  to  control  system  compensation,  while  avoiding  the 
painstaking  labor  of  trial-and-error  analysis  inherent  in 
that  technique.  Several,  practical  engineering  examples 
have  been  presented  to  demonstrate  the  superiority  of  the 
graphical  technique.  To  date,  no  other  package  is  known  to 
offer  the  fully  interactive  and  comprehensive  capabilities 
of  the  parameter  plane  program. 
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9 
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By  itself  the  program  allows  one  to  design  a  control 
system  compensation  model  for  most  systems.  However,  for 
some  lightly  damped  systems  containing  mechanical  resonances, 
the  amount  by  which  zeta  or  omega  are  incremented  in  the 
parameter  plane  equations  may  be  so  large  as  to  not  detect 


the  resonance  peaks.  This  information  would  be  available  from 
either  a  root-locus  or  Bode  analysis.  For  still  other 
systems,  one  might  be  interested  in  the  way  the  roots  of  the 
characteristic  equation  extend  from  the  open  loop  poles  and/or 
zeros.  Since  the  parameter  plane  equations  are  calculated 
using  only  the  characteristic  equation,  no  knowledge  of  open 
loop  poles  or  zeros  is  available.  Again,  a  root-locus  method 
would  reveal  this  information.  Incorporated  into  one 
comprehensive  package  which  includes  Bode  and  root-locus 
analyses,  the  program  provides  the  capability  to  investigate 
the  entire  gamut  of  linear  control  system  architecture. 

A  basis  for  further  investigation  involves  plotting  the 
parameter  plane  contours  for  systems  that  are  non-linear 
in  the  alpha  and  beta  terms--i.e.,  those  systems  which  contain 
alpha-beta  product  terms.  Although  the  recursion  technique 
used  in  this  text  has  distinct  advantages  over  the  matrix 
approach  for  the  linear  case,  as  pointed  out  earlier,  the 
matrix  technique  would  be  the  method  of  choice  for  the  non¬ 
linear  case. 
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APPENDIX  £ 

PARAMETER  PLANE  PROGRAM 


C 

C 


C 

C 

100 


101 

C 

C 

102 


C 

C 

103 


104 

105 
C 

C 

106 


C 

C 

107 


C 


SUBROUTINE  LPARAM 
DIMENSION  AC 350 ) ,  BC350), 

X  AGC  350 ) ,  BGC 350 ) ,  AJ(IOO),  BJ(IOO),  CJC100), 

X  ZETA(IOO),  SIGMA(IOO),  W(100),  ZW(IOO) 

INITIAL  ASSIGNMENTS 

CHARACTERX4  SCHAR/'S  '/,YES/'Y  '/,NOO/'N  */, BLANK/'  '/, 
X  ENCHAR/'E  •/, WNCHAR/ ' WN '/, NDCHAR/ ' ND'/ , NOCHAR/ ' NO '/ , 
X  AJCHAR/ ' A J '/ , BJCHAR/ ' B J  V , CJCHAR/ ' CJ '/ , NSCHAR/ ' NS ’/ , 

X  NZCHAR/ • NZ • / , ZWCHAR/ • ZW  * / , NWCHAR/ • NW • / , PRCHAR/ *  PR ' / , 

X  NCCHAR/ ' NC'/ ,  TICHAR/  '  TI '/ 

CHARACTERX4  TABLE,  GRAPH,  CHANGE,  REPLY,  OPT,  LABEL (9) 
COMMON  /SAVE/  LABEL,  WN,  ND,  N02,  NC,  CJ,  AJ,  BJ, 

X  NZ,  ZETA,  NS,  SIGMA,  NW,  W,  NZW,  ZW, 

X  XMIN,  XMAX,  YMIN,  YMAX 

DATA  ENTRY  FROM  FILE  OR  CONSOLE? 

MINMAX  =  1 
GRD  =  0. 

CHANGE  =  BLANK 
CALL  EXCMSC 'CLRSCRN' ) 

WRITEC6 , 500 ) 

CALL  REA DC  (REPLY) 

IF  (REPLY  .NE.  'D')  GOTO  101 
CALL  GETIT 
MINMAX  =  0 
GO  TO  200 
CONTINUE 


GET  CURVE  TITLE 

CONTINUE 

CALL  EXCMS( 'CLRSCRN' ) 

WRITE( 6 , 501 ) 

CALL  READL  (LABEL) 

CALL  ASTER  ( LABEL , LABEL ) 

IF  (  CHANGE  .EQ.  TICHAR  )  GO  TO  200 

GET  STARTING  VALUE  OF  WN 

CONTINUE 
WRITE(6 , 502) 

CALL  READR  (WN) 

IF  (WN)  104,104,105 
WRITE(6 , 503) 

GO  TO  103 

IF  (  CHANGE  .EQ.  WNCHAR)  GOTO  200 

GET  THE  NUMBER  OF  DECADES  CONSIDERED 

CONTINUE 
WRITE( 6 , 504 ) 

CALL  READI  (ND) 

IF  (  CHANGE  .EQ.  NDCHAR  )  GOTO  200 

GET  THE  ORDER  OF  THE  CHARACTERISTIC  EQN 

CONTINUE 
WRITE( 6,505) 

CALL  READI  (N02) 

NC  =  N02+1 

IF  (CHANGE  .EQ.  NOCHAR  )  GOTO  200 


GET  THE  NUMBER  OF  CONSTANT  ZETA  CURVES 


i 


I 


§ 


I 


C 

C 


108 


C 

C 


109 


110 


C 

C 


111 


C 

C 


112 


113 


C 

C 


114 


C 

C 


115 


116 


C 

C 


117 


C 

C 


118 


119 


C 

c 

c 

c 


CONTINUE 

CALL  EXCMS  ('CLRSCRN') 
WRITE(6 ,  506 ) 

CALL  READI  (NZ) 

IF  (NZ  .LT.  1)  GOTO  110 


GET  THE  VALUES  OF  ZETA 

WRITE(6 , 507 ) 

DO  110  I  =  1 , NZ 
WRITE(6, 508)1 
CALL  READR  (ZETA(I)) 

IF  (ZETA(I)  .LT.  0.  .OR.  ZETA(I)  .GT.  1 
IF  (ZETA(I)  .LT.  0.  .OR.  ZETA(I)  .GT.  1 
CONTINUE 

IF  (  CHANGE. EQ.  NZCHAR  )  GOTO  200 


WRITE(6 , 509 ) 
GO  TO  109 


CONTINUE 
CALL  EXCMS  ( 'CLRSCRN' ) 
WRITEC6 , 510 ) 

CALL  READI  (NS) 

IF  (NS  .LT.  1)  GOTO  113 


GET  THE  NUMBER  OF  CONSTANT  SIGMA  CURVES 


GET  THE  VALUES  OF  SIGMA 

DO  113  I  =  1 , NS 
WRITE(6 ,  511 )  I 
CALL  READR  (SIGMA(I)) 

IF  (SIGMA(I)  .LT.  0.)  WRITE  (6,512) 

IF  (SIGMA(I)  .LT.  0.)  GO  TO  112 
CONTINUE 

IF  (  CHANGE  .EQ.  NSCHAR  )  GOTO  200 


CONTINUE 
CALL  EXCMS  ( 'CLRSCRN' ) 
WRITE(6,513) 

CALL  READI  (NW) 

IF  (NW  .LT.  1)  GOTO  116 


GET  THE  NUMBER  OF  CONSTANT  WN  CURVES 


GET  THE  WN  VALUES 

WNMAX  =  WN*10*XND 
DO  116  I  =  1 , NW 
WRITE(6,514)  I 
CALL  READR  (W(I)) 

IF  (W(I)  .LT.  WN  .OR.  W(I)  .GT.  WNMAX)  WRITE  (6,515)  WN, WNMAX 
IF  (W(I)  .LT.  WN  .OR.  W(I)  .GT.  WNMAX)  GO  TO  115 
CONTINUE 

IF  (  CHANGE  .EQ.  NWCHAR  )  GOTO  200 


GET  THE  NUMBER  OF  CONSTANT  ZETAXWN  CURVES 

CONTINUE 

CALL  EXCMS  ('CLRSCRN') 

WRITE( 6 , 516 ) 

CALL  READI  (NZW) 

IF  (NZW  .LT.  1)  GOTO  119 


GET  THE  Z*WN  VALUES 

DO  119  I  =  1 , NZW 
WRITE(6 , 517  )  I 
CALL  READR  (ZW(I)) 

IF  (ZW(I)  .LE.  0.)  WRITE  (6,518) 

IF  (ZW(I)  .LE.  0.)  GO  TO  118 
CONTINUE 

IF  (  CHANGE  .EQ.  ZWCHAR  )  GOTO  200 
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c 

c 

120 


121 


C 

C 

122 


123 


C 

C 

124 


125 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


GET  CONSTANT  COEFFICIENTS 

CONTINUE 

CALL  EXCMS  ('CLRSCRN') 

WRITE(6 ,519) 

DO  121  I  =  NC,1,-1 
II  =  1-1 

WRITE! 6 , 520 )  II, I 
CALL  READR  CCJ(I)) 

CONTINUE 

IF  (  CHANGE  .EQ.  CJCHAR  )  GOTO  200 

GET  ALPHA  COEFFICIENTS 

CONTINUE 

CALL  EXCMS  ('CLRSCRN') 

WRITEC6 , 521 ) 

DO  123  I  =  NC,1,-1 
II  =  1-1 

WRITE(6,522)  II, I 
CALL  READR  (AJCI)) 

CONTINUE 

IF  (CHANGE  .EQ.  AJCHAR  )  GOTO  200 

GET  BETA  COEFFICIENTS 

CONTINUE 

CALL  EXCMS  ('CLRSCRN') 

WRITE(6 , 523) 

DO  125  I  =  NC,1,-1 
II  =  1-1 

WRITE(6 , 524)  11,1 
CALL  READR  (BJ(I)) 

CONTINUE 
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C  REVIEW  ENTRIES 

200  CONTINUE 

CALL  EXCMS  C'CLRSCRN') 

WRITEC6 , 525) 

CALL  REA DC  (REPLY) 

IF  (REPLY  .NE.  YES)  GOTO  209 
CALL  EXCMS  ('CLRSCRN') 

WRITE(6,526) 

WRITE(6  > 527 )  (LABEL(I),  1  =  1,  9) 
WRITE(6,528) 

WRITE(6 , 529 )  NO, N02, NZ, NS, NW, NZW 
WRITE(6 , 530 )  WN 
WRITE(6 , 531 ) 

IF  (NZ  .LT.  1)  GOTO  201 

MRITE(6 , 532)  (ZETA(M) ,M=1 , NZ) 

GOTO  202 

201  WRITE(6 , 533) 

202  CONTINUE 
WRITE(6,534) 

IF  (NS  .LT.  1)  GoTO  203 

WRITE(6,532)  (SIGMA(M),M=1,NS) 

GOTO  204 

203  WRITE(6,533) 

204  CONTINUE 
WRITE(6 , 535) 

IF  (NW  .LT.  1)  GOTO  205 
WRITE(6,532)  (W(M) ,M=1 , NW) 

GOTO  206 

205  WRITE(6,533) 

206  CONTINUE 
WRITE(6 , 536 ) 

IF  (NZW  .LT.  1)  GOTO  207 

WRITE(6 , 532)  (ZW(M),M=1,NZW) 

GOTO  208 

207  WRITE(6,533) 

208  CONTINUE 
WRITE  (6,576) 

CALL  EXCMS ( 'CLRSCRN*) 

WRITE(6 , 537 ) 

WRITE(6,532)  (CJ(N) , N=NC, 1 , -1 ) 
WRITE(6,538) 

WRITE(6 , 532)  ( AJ(N) , N=NC, 1 , -1 ) 
WRITE(6,539) 

WRITE(6 , 532)  ( BJ( N) , N=NC, 1 , -1 ) 
WRITE(6,540) 

WRITE(6 , 541 )  XMIN,  XMAX,  YMIN,  YMAX 

209  CONTINUE 
WRITE(6 , 542) 

CALL  REA DC  (CHANGE) 

IF  (CHANGE  .NE.  YES)  GOTO  210 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 

c 

c 

c 
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CHANGE  ROUTINE 


CALL  EXCMS  CCLRSCRN') 
WRITE(6,543) 

CALL  READC  (CHANGE) 

IF  (CHANGE  .EQ.  TICHAR) 
IF  (CHANGE  .EQ.  WNCHAR) 
IF  (CHANGE  .EQ.  NDCHAR) 
IF  (CHANGE  .EQ.  NOCHAR) 
IF  (CHANGE  .EQ.  NZCHAR) 
IF  (CHANGE  .EQ.  NSCHAR) 
IF  (CHANGE  .EQ.  NWCHAR) 
IF  (CHANGE  .EQ.  ZWCHAR) 
IF  (CHANGE  .EQ.  CJCHAR) 
IF  (CHANGE  .EQ.  AJCHAR) 
IF  (CHANGE  .EQ.  BJCHAR) 
IF  (CHANGE  .EQ.  ENCHAR) 
IF  (CHANGE  .EQ.  PRCHAR) 
IF  (CHANGE  .EQ.  NCCHAR) 
GOTO  209 
CONTINUE 
WRITE(6 , 544) 

CALL  READC  (TABLE) 
WRITE(6 , 545) 

CALL  READC  (GRAPH) 


IF  (NZ)  213,213,212 
G  =  (10.xxND)xx(l./299.) 
CONTINUE 
Z  =  1.0E-60 


GOTO  102 
GOTO  103 
GOTO  106 
GOTO  107 
GOTO  108 
GOTO  111 
GOTO  114 
GOTO  117 
GOTO  120 
GOTO  122 
GOTO  124 
GOTO  405 
GOTO  100 
GOTO  210 
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CONSTANT  ZETA  PLOTS 


300 

C 

C 


C 


C 


301 

302 


303 
C 

304 


C 


305 

306 
C 

C 


IF  (NZ)  309,309,300 
CALL  EXCMS  ( 'CLRSCRN' ) 

IF  (TABLE  .EQ.  YES  .OR.  GRAPH  .EQ.  YES)  WRITE  (6,546) 


DO  308  M=1 , NZ 

IF  (TABLE  .EQ.  YES)  WRITE  (6,547) 

J  =  0 
R  =  0. 

WNA  =  WN 

DO  306  L-1,300 
D1  =  0.0 
D2  =  0.0 
Cl  =  0.0 
C2  =  0.0 
B1  =  0.0 
B2  =  0.0 

DO  303  N=1 , NC 
K  -  N~1 

IF  (K)  302,301,302 
U  =  0.0 
U1  =  -1.0 

U2  =  2 . 0*ZETA(M)*U-U1 
D1  =  (-1 . 0)**K*CJ(N)*WNAXXKXU1+D1 
D2  =  (-1 .0)XXKXCJ(N)XWNAXXKXU+D2 
Cl  =  (-1 .0)XXKXBJ(N)XWNAXXKXU1+C1 
C2  =  (-1 .0)XXKXBJ(N)XWNAX*KXU+C2 
B1  =  (-1 .0)XX«XAJ(N)XWNAXXKXU1+B1 
B2  =  (-1.0)X*KXAJ(N)XWNAXXKXU+B2 
U1  =  U 
U  =  U2 
CONTINUE 

IF  (ABS(B1XC2-B2XC1)-Z)  306,306,304 
J  =  J  +  l 
R  =  R+l . 

A( J )  =  (C1XD2-C2XD1)/(B1XC2-B2XC1) 

B( J )  =  (B2XD1-B1XD2)/(B1XC2-B2XC1) 

IF  (TABLE  .NE.  YES)  GO  TO  306 

WRITE(6 , 548 )  A(J),  B(J),  WNA,  ZETA(M) 

IF  (R/10.  -  J/10)  306,305,306 

CALL  ROOTS  (A(J),  B(J),  AJ,  BJ,  CJ ,  N02) 

CALL  EXCMS( 'CLRSCRN') 

WRITE  (6,547) 

WNA  =  GXWNA 

CALL  EXCMS( 'CLRSCRN* ) 

IF  (J  .GT.  0)  GOTO  307 
WRITE(6 , 549 ) 

GOTO  308 

IF  (GRAPH  .EQ.  YES)  CALL  PLOTD( A , B, J , . FALSE . , 

LABEL,  ' ALPHAS ' , ' BETA$  ' , 
MINMAX,'  Z=$ ' ,ZETA(M) , 
XMI N , XMAX , YMI N , YMAX , GRD ) 

IF  (GRAPH  .EQ.  YES)  GRD  =  GRD+1 . 

CONTINUE 

CONTINUE 
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C 

c 

325 


C 

C 


C 


326 
C 

327 

328 

329 
C 

330 


331 


332 
C 

333 


CONSTANT  SIGMA  PLOTS 

IF  (NS)  335,335,325 
CALL  EXCMSC 'CLRSCRN' ) 

IF  (TABLE  .EQ.  YES  .OR.  GRAPH  . EQ.  YES)  WRITE  (6,550) 

XINC  =  (XMAX  -  XMIN)/299. 

YIN C  =  (YMAX  -  YMIN)/299 . 

DO  334  M=1 , NS 
XPT  =  XMIN 
YPT  =  YMIN 
DD  =  CJ(1) 

CC  =  BJ( 1 ) 

BB  =  A J ( 1 ) 

J  =  0 
R  =  0.  ■ 

DO  326  N=2, NC 
K  =  N-l 

DUMMY5=SIGMA(M)**K 
DD  =  (-1 . 0)**KXCJ(N)*DUMMY5+DD 
CC  =  (-1.0)**K*BJ(N)*DUMMY5 +CC 
BB  =  (-1.0)X*KXAJ(N)*DUMMY5+BB 
CONTINUE 

IF  (CC  .EQ.  0.  .AND.  BB  .EQ.  0.)  GOTO  334 
IF  (CC)  327,327,330 
DO  329  L=l,300 
J  =  J  +  l 
R  =  R+l. 

A( J )  =  XPT 

B(J)  =  (-BB*A(J)-DD)/CC 
IF  (TABLE  .NE.  YES)  GOTO  329 
WRITE(6 , 552)  A(J),  B(J),  SIGMA(M) 

IF  (R/10.  -  J/10)  329,328,329 

CALL  ROOTS  (A(J),  B(J),  AJ,.  BJ,  CJ,  N02) 

CALL  EXCMS  ('CLRSCRN') 

WRITE(6,550) 

XPT  =  XPT  +  XINC 
GO  TO  333 

DO  332  L=l,300 
J  =  J+I 
R  =  R+l. 

B( J )  =  YPT 

A( J )  =  (-CC*B(J)-DD)/BB 
IF  (TABLE  .NE.  YES)  GOTO  332 
WRITE(6 , 552 )  A(J),  B(J),  SIGMA (M) 

IF  (R/10.  -  J/10)  332,331,332 

CALL  ROOTS  (A(J),  B(J),  AJ,  BJ,  CJ ,  N02) 

CALL  EXCMS  ('CLRSCRN') 

WRITE(6 , 550 ) 

YPT  =  YPT  +  XINC 

CALL  EXCMS( 'CLRSCRN') 

IF  (GRAPH  .EQ.  YES)  CALL  PLOTD( A , B, J , . FALSE . , 

LABEL,  'ALPHA$' , 'BETAS' , 
MINMAX, '  S=$ ' , SIGMA(M) , 
XMIN, XMAX, YMIN, YMAX, GRD) 

IF  (GRAPH  .EQ.  YES)  GRD  =  GRD+I . 

CONTINUE 

CONTINUE 
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350 


C 

C 

C 


c 


c 


351 

352 


353 
C 

354 


355 


356 

C 

C 


C 

357 

X 

X 

X 

C 

C 

358 
C 

359 
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CONSTANT  ZETA-OMEGA  PLOTS 
IF  (NZW)  359,359,350 
CALL  EXCMSC 'CLRSCRN' ) 

IF  (TABLE  . EQ .  YES  .OR.  GRAPH  .EQ.  YES)  WRITE  (6,553) 

XWN  =  WN 

DO  358  M=1,NZW 

IF  (TABLE  .EQ.  YES)  WRITE  (6,554) 

J  =  0 
R  =  0. 

AZETA  »  1./300. 

DO  356  L^l, 300 

XWN  =  ZW(M)/AZETA 

D1  a  0.0 

D2  =  0.0 

Cl  =  0.0 

C2  =  0.0 

B1  =  0.0 

B2  a  0.0 

DO  353  N=1,NC 

X  a  N-l 

IF  (K)  352,351,352 
Q1  =  0.0 
Q  =  -1.0/XWNKX2 
D2  =  CJ(N)XQ1+D2 
C2  =  BJ ( N )*Q1+C2 
B2  =  AJ(N)XQ1+B2 
D1  =  CJ(N)XQ+D1 
Cl  =  BJ(N)*Q+C1 
B1  =  AJ(N)*Q+B1 
Q2  =  -2 . 0*ZW(M)xQl-XWNX*2*Q 
Q  =  Q1 
Q1  =  Q2 

IF  (ABS(B1*C2-B2XC1)-Z)  35(5,356,354 
J  =  J+l 
R  =  R+l. 

A(  J )  =  (C1XD2-C2XD1 )/( B1*C2-B2*C1 ) 

B(J)  =  (B2XD1-B1XD2)/(B1XC2-B2*C1) 

IF  (TABLE  .NE.  YES)  GO  TO  356 
WRITE(6 , 552)  A(J),  B(J),  ZW(M) 

IF  (R/10.  -  J/10)  356,355,356 

CALL  ROOTS  (A(J),  B(J),  AJ,  BJ,  CJ,  N02) 

CALL  EXCMS( 'CLRSCRN*) 

WRITE  (6,554) 

AZETA  =  AZETA+U  ./300. ) 

CALL  EXCMS( 'CLRSCRN') 

IF  (J  .GT.  0)  GOTO  357 
WRITE(6 , 549 ) 

GOTO  358 

IF  (GRAPH  .EQ.  YES)  CALL  PLOTD( A, B, J , . FALSE . , 

LABEL,  'ALPHAS ' , ' BETAS  * , 
MINMAX, ' ZW=S ' , ZW( M) , 
XMIN, XMAX,YMIN, YMAX.GRD) 

IF  (GRAPH  .EQ.  YES)  GRD  =  GRD+1. 

CONTINUE 

CONTINUE 
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CONSTANT  OMEGA  PLOTS 

IF  (NW)  385,385,376 
CALL  EXCMSC 'CLRSCRN' ) 

IF  (TABLE  .EQ.  YES  .OR.  GRAPH  .EQ.  YES)  WRITE  (6,555) 


DO  384  M=I , NW 

IF  (TABLE  . EQ .  YES)  WRITE  (6,547) 

J  =  0 
R  =  0. 

AZETA  =0.0 

DO  382  L=1 , 300 
D1  =  0.0 
D2  =  0.0 
Cl  =  0.0 
C2  =  0.0 
B1  =  0.0 
B2  =  0.0 

DO  379  N=1 , NC 
K  =  N-l 

IF  (K)  378,377,378 
U  =  0.0 
U1  =  -1.0 

U2  =  2.0XAZETAXU-U1 
D1  =  (-1.0)XXKXCJ(N)XW(M)XKKXU1+D1 
D2  =  (-1.0)X*K*CJ(N)XW(M)XXKXU+D2 
Cl  =  (-1.0)X*KXBJ(N)XW(M)XXKXU1+C1 
C2  =  (-1 .0)XXKXBJ(N)XW(M)XXK»U+C2 
B1  =  (-1.0)XXKXAJ(N)XW(M)XXKXU1+B1 
B2  =  ( -1 . 0 )X*KXAJ (N )XW(M) XXKXU+B2 
U1  =  U 
U  =  U2 
CONTINUE 

IF  (ABSCB1XC2-B2XCD-Z)  382,382,380 
J  =  J  +  l 
R  =  R+l. 

A ( J )  =  (ClxD2-C2»Dl )/( B1XC2-B2XC1 ) 

B( J )  =  (B2XD1-B1XD2)/(B1XC2-B2XC1) 

IF  (TABLE  .NE.  YES)  GO  TO  382 
WRITE(6 , 548 )  A(J),  B(J),  W(M),  AZETA 
IF  (R/10.  -  J/10 )  382,381,382 
CALL  ROOTS  (A(J),  B(J),  AJ,  BJ,  CJ,  N02) 

CALL  EXCMS( 'CLRSCRN' ) 

WRITE  (6,547) 

AZETA  =  AZETA+U./299.) 

CALL  EXCMS( 'CLRSCRN' ) 

IF  (J  .GT.  0)  GOTO  383 
WRITE(6 , 549 ) 

GOTO  384 

IF  (GRAPH  .EQ.  YES)  CALL  PLOTD(A, B, J, . FALSE. , 

LABEL,  '  ALPHA$ ' , ' BETAS ' , 
MINMAX,'  W=$ ' , W(M) , 
XMIN,XMAX,YMIN,YMAX,GRD) 


CONTINUE 

CONTINUE 


X 

X 

X 


IF  (GRAPH  .EQ. 


400 


401 


IF  (GRAPH  .EQ. 

GRD  =  0. 

CALL  EXCMS( 'CLRSCRN' ) 

IF  (OPT  .EQ.  YES)  CALL  DONEPL 
YES)  STOP 


YES)  CALL  PLOTD(0. ,0. ,0. , .TRUE. , 
’  $','  $', 
0.,'  $',-9.7531, 
0.,0.,0.,0.,GRD) 
YES)  GRD  =  GRD+1 . 


IF  (OPT  .EQ. 

WRITE(6 , 557 ) 
WRITE(6,556 ) 

WRITEC6 , 557 ) 

WRITE( 6 , 558 ) 

WRITE(6 , 559 ) 

WRITE(6 , 560 ) 
WRITE(6, 561 ) 

WRITE(6 , 562) 

WRITE(6 , 563 ) 

WRITE(6 , 557 ) 

CALL  READI  (IANS) 

IF  (IANS  .GT.  6 
GO  TO  (100,  404, 


OR.  IANS  .LT.  1)  GOTO  400 
401,  402,  403,  405)  IANS 


ROOT  FINDER  OPTION 

CALL  EXCMS( 'CLRSCRN') 

WRITE(6 , 564) 

WRITE(6 , 565) 

CALL  READR  (ALPHA) 

WRITE(6 , 566 ) 

CALL  READR  (BETA) 

CALL  ROOTS  (ALPHA,  BETA,  AJ,  BJ,  CJ,  N02) 


GO  TO  400 


402 


CALL  SAVIT 
GO  TO  400 


SAVE  OPTION 


403 


CREATE  DISSPLA  METAFILE  OPTION 

CALL  EXCMS( 'CLRSCRN') 

WRITE(6 , 567 ) 

WRITE(6 , 568 ) 

WRITE(6 , 569 ) 

CALL  READC  (OPT) 

IF  (OPT  .NE.  YES)  GOTO  400 
CALL  EXCMS( 'CLRSCRN') 

WRITE(6,570) 

CALL  DONEPL 


CALL  META 
GO  TO  211 


COMPRS  =  SUBROUTINE  TO  LET  DONEPL  FINISH 


SAME  PROBLEM  OPTION 

404  WRITE(6 , 571 ) 

CALL  READI  (MINMAX) 

IF  (MINMAX  .EQ.  1)  GO  TO  200 
CALL  EXCMS( 'CLRSCRN' ) 

WRITE(6 ,572) 

CALL  READR  (XMIN) 

WRITE( 6 ,573) 

CALL  READR  (XMAX) 

WRITEC  6 , 574 ) 

CALL  READR  (YMIN) 

WRITEC  6 , 57  5 ) 

CALL  READR  ( YMAX) 

GO  TO  200 

405  CONTINUE 
RETURN 


oooooooon  oooooooooo 


wv 


FORMATS  STATEMENTS 

F0RMATC5X, 'THIS  IS  THE  INTERACTIVE  PARAMETER  PLANE  PROGRAM...*,/, 

X  5X, 'THE  USER  WILL  BE  PROMPTED  FOR  VARIOUS  INPUTS.’,////, 

X  5X, 'WILL  YOU  BE  ENTERING  DATA  FROM  A  CONSOLE  OR  DATAFILE?', 

X  / ,  15X,  ' "D"  OR  "C") 

FORMAT!///, '  ENTER  TITLE  TO  APPEAR  FOR  THIS  FAMILY  OF  CURVES') 
FORMAT!/, IX, 'WHAT  IS  THE  STARTING  VALUE  OF  OMEGAN  !WN>0.0)?’) 
FORMAT!/, IX, ’WN  MUST  BE  GREATER  THAN  ZERO  -  TRY  AGAIN') 

FORMAT!/, IX, ’HOW  MANY  DECADES  PAST  WN  ARE  DESIRED?  ! ND) ’ ) 
FORMAT!/,’  WHAT  IS  THE  ORDER  OF  THE  CHARACTERISTIC  EQUATION?!NO) ’ ) 
FORMAT!/, IX, ’HOW  MANY  CONSTANT  ZETA  CURVES  ARE  DESIRED?  !NZ)’) 
FORMAT!/,'  ENTER  THE  VALUES  OF  ZETA  TO  BE  USED  IN  COMPUTATION...’) 
FORMAT ! 5X, 'ZETA! ' , 12, ’ )  =  ?•) 

F0RMAT!5X, 'ZETA  MUST  LIE  BETWEEN  0  AND  1,  INCLUSIVE  -  TRY  AGAIN’) 
FORMAT!/,'  HOW  MANY  CONSTANT  SIGMA  !REAL  ROOT)  CURVES  ARE  DESIRED? 
X  !NS) ' ) 

F0RMATL5X, ’SIGMA! *, 12, ')  =  ?') 

FORMAT! 5X, ’NEGATIVE  SIGMA  MEANS  POSITIVE  REAL  ROOT  -  TRY  ANOTHER’) 
FORMAT!/, IX, ’HOW  MANY  CONSTANT  WN  CURVES  ARE  DESIRED?  CNW)') 
F0RMAT!5X, ’W! • ,12, ')  =  ?’) 

FORMAT  ! 5X, ’ WN  NOT  WITHIN  PLOTTABLE  RANGE', 

X  / , 5X, ’YOUR  USABLE  RANGE  IS' 

X  / , IOX, FI 0 . 2, '  TO  ' , F10 . 2) 

FORMAT!/, IX, 'HOW  MANY  CONSTANT  Z*WN  CURVES  ARE  DESIRED?  INZW)') 
F0RMAT!5X, 'ZW! ',12, ')  =  ?') 

F0RMATI5X, ’NON-POSITIVE  Z-WN  MEANS  POSITIVE  ROOT  -  TRY  ANOTHER') 
FORMAT!/, IX, 'ENTER  THE  CONSTANT  COEFFICIENTS...’) 

FORMAT ! 5X, ’ - S*x ’ , 12, ’ -  CJ!’,I2,’)  =  ?') 

FORMAT!/, IX, ’ENTER  THE  ALPHA  COEFFICIENTS . . . ’ ) 

FORMAT ! 5X, ’ ---SXX ’ , 12, ' -  AJ!’,I2,’)  =  ?’) 

FORMAT!/, IX, ’ENTER  THE  BETA  COEFFICIENTS...’) 

FORMAT ! 5X,  ’ - SXX’,12,’—  BJ!’,I2,’)  =  ?’) 

FORMAT!/,'  WANT  TO  REVIEW  YOUR  ENTRIES  BEFORE  RUNNING?  !Y/N)') 
FORMAT  !/, 10X, ’  GRAPH  TITLE’) 

FORMAT  ! IX, 9A4 ) 

FORMAT  ! /, 8X, 2HND, 8X, 2HN0, 8X, 2HNZ, 8X, 2HNS, 8X, 2HNW, 7X, 3HNZW) 

FORMAT  I6I10) 

FORMAT  !/,10X, ’INITIAL  VALUE  OF  OMEGA  =  ’,F10.5) 

FORMAT  !/, 10X, '  ZETA  ’) 

FORMAT  I8E10.3) 

FORMAT  !  IX,  ' . NO  VALUE . ’  ) 

FORMAT  !/,10X, 'SIGMA  ') 

FORMAT  !/, 10X, '  W  ') 

FORMAT  ! /, 1  OX, 'ZW  ' ) 

FORMAT  I/.10X, 'CONSTANT  COEFFICIENTS  IN  DECENDING  ORDER') 

FORMAT  !/,10X, ’ALPHA  COEFFICIENTS  IN  DECENDING  ORDER’) 

FORMAT  ! / , 10X, ’BETA  COEFFICIENTS  IN  DECENDING  ORDER’) 

FORMAT  (/, IOX, ’XMIN  XMAX  YMIN  YMAX ’ ) 

FORMAT  ! IX, 4E10 . 3 ) 

FORMAT!/, ’  WANT  TO  MAKE  ANY  CHANGES?  !Y/N)’) 

FORMAT!/,’  WHAT  VARIABLE/AREA  DO  YOU  WISH  TO  CHANGE?’,//, 

X5X. ’TITLE . TI  OMEGA  START  ..  WN  #  DECADES _ ND’,/ 


...  V..'  . 


V 


X5X, 'ORDER . NO  CONST  ZETA...NZ  CONST  SIGMA.. NS',/ 

X5X, 'CONST  WN . NW  CONST  Z-WN . . . ZW' ,// 

X5X, 'CONST  TERMS.. CJ  ALPHA  TERMS .. AJ  BETA  TERMS ... BJ ' 

X5X, ' END . E  NEW  PROBLEM. .PR  NO  CHANGE _ NC’,///, 

XIX, 'ENTER  TWO  DIGIT  CODE...') 

FORMATC/,'  DO  YOU  WANT  TABULATED  DATA  ON  THE  SCREEN?  CY/N)') 
FORMATC/,'  DO  YOU  WANT  THE  GRAPHS  ON  THE  TERMINAL?  (Y/N)') 

FORMAT  (1H1,10X, 'CONSTANT  ZETA  CURVES') 

FORMATC/, 10X, 'ALPHA  BETA  OMEGA  ZET 

FORMAT  (4E16.5) 

FORMATC/,'  DUE  TO  PLOT  RESTRICTIONS,  A  COMPLETE  GRAPH  CANNOT  BE 
XTPUT . ' ) 

FORMAT  C1H1,10X, 'CONSTANT  SIGMA  CURVES') 

FORMAT  C/,10X, 'ALPHA  BETA  SIGMA') 

FORMAT  C3E16.5) 

FORMAT  C 1H1 , 10X, 'CONSTANT  ZETA-OMEGA  CURVES’) 

FORMAT  C/,10X, 'ALPHA  BETA  ZETA-OMEGA') 


OMEGA 


FORMAT  C 1H1 , 10X, 'CONSTANT  ZETA-OMEGA  CUR\ 
FORMAT  C/,10X, 'ALPHA  BETA 

FORMAT  C1H1.10X, 'CONSTANT  OMEGA  CURVES’) 


ZETA-OMEGA' ) 


FORMATC 5X, 
FORMAT C5X, 
FORMAT  C  5X, 
FORMAT C5X, 
FORMATC 5X, 
FORMATC 5X, 
FORMAT C5X, 
FORMAT C5X, 
FORMAT  C  5X, 

Y  CY 

FORMATC/// 


FORMATC/, 5X. 'ENTER  THE  BETA  VALUE 


FORMAT  C  5X, 

X  5X, 

X  5X, 

X  5X, 

X  5X, 

FORMATC 5X, 
X  5X , 

Y  cv 

FORMATC/// 
X  15X 


I  OPTION  NO. I 


OPTION 


111  NEW  PROBLEM 

I  2  I  SAME  PROBLEM 

I  3  |  ROOT  FINDER 

I  4  |  SAVE  DATA 

I  5  |  SAVE  GRAPH  IN  DISSPLA  METAFILE 

I  6  I  RETURN  TO  MAIN  MENU 

ENTER  AN  ALPHA-BETA  PAIR,  AND  THE  ROOTS  OF  YOUR 
SYSTEMS  CHARACTERISTIC  EQUATION  WILL  BE  RETURNED 
5X, 'ENTER  THE  ALPHA  VALUE 


YOU  NOW  HAVE  THE  OPTION  OF  STORING  THE  LAST  SET  OF  '/, 
CURVES  IN  A  DISSPLA  METAFILE.  THIS  ALLOWS  RETRIEVAL’/, 


OF  DATA  AT  A  LATER  TIME  FOR  ROUTING  TO  ANY  OF 
SEVERAL  OUTPUT  DEVICES  CTEK618,  3800  LASER  PRINTER, 
VERSATEC  PLOTTER,  ETC.) 

IF  YOU  CHOOSE  THIS  OPTION,  THE  PROGRAM  MUST  BE 
TERMINATED  -  THIS  CANNOT  BE  AVOIDED  WITHOUT 
CATASTROPHIC  RESULTS.  ' 

5X, 'DO  YOU  WISH  TO  USE  THIS  OPTION? 

•"Y"  OR  "N" 


FORMATC/////, 5X, 'IF  YOU  WISH  GRAPHIC  OUTPUT,  TYPE: 

X  15X, '"DISSPOP" 

X  5X , ' AND  FOLLOW  THE  INSTRUCTIONS  .  .  . 

FORMAT  C///,'  AUTOSCALE  OR  USER-DEFINED  LIMITS  FOR  CURVES?’, 

X/,'  1 =AUTOSCAL  E ;  0  =  USER-DEFINED' ) 

FORMAT  C//,'  INPUT  MINIMUM  VALUE  FOR  X  CX-MIN)') 

FORMAT  C/.'  INPUT  MAXIMUM  VALUE  FOR  X  CX-MAX)') 

FORMAT  C/,'  INPUT  MINIMUM  VALUE  FOR  Y  CY-MIN)') 

FORMAT  C/,'  INPUT  MAXIMUM  VALUE  FOR  Y  CY-MAX)') 

FORMAT  C////////////////) 

END 
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SUBROUTINE  PLOTD  —  GRAPHS  WILL  BE  PRODUCED  ON  UPRIGHT  11  X  14 
PAGE  WITH  9  INCH  AXES.  IF  USER  SELECTS  'AUTOSCALE'  FEATURE, 
SUBROUTINE  PLD010  (INTERNAL  TO  PLOTD)  FINDS  MIN  AND  MAX  FOR  EACH  = 
AXIS  AND  SCALES  ACCORDINGLY.  FORMAT. 

CALL  PLOTDCXDATA,  YDATA,  NNPTS,  EJECT,  LABEL,  XL ABEL ,  YLABEL , 

X  MINMAX,  CRVTTL,  CRVNUM,  XMIN,  XMAX,  YMIN,  YMAX) 

WHERE. 

XDATA  IS  A  REALX4  ARRAY  DIMENSIONED  AT  LEAST  | NNPTS! 
CONTAINING  THE  X  ORDINATE  VALUES, 

YDATA  IS  A  REALX4  ARRAY  DIMENSIONED  AT  LEAST  INNPTSI 
CONTAINING  THE  Y  ORDINATE  VALUES, 

NNPTS  IS  AN  INTEGERX4  SCALAR  DESIGNATING  THE  NUMBER  OF  = 
POINTS  TO  BE  PLOTTED.  THE  NUMBER  OF  POINTS  IS 
ABS(NNPTS) .  NNPTSCO  MEANS  PLOT  POINTS  ONLY. 

EJECT  IS  A  L0GICALX4  VARIABLE  OR  CONSTANT  INDICATING  = 
WHETHER  A  PAGE  EJECT  IS  REQUIRED  FOLLOWING  THE  = 
CURRENT  CURVE.  THIS  ALLOWS  MULTIPLE  CURVES  ON  ONE  = 
SET  OF  EXES.  PAGE  EJECT  WILL  OCCUR  FOR  NEXT  GRAPH  = 
AFTER  EJECT  HAS  BEEN  SET  TO  .TRUE. 

LABEL  IS  A  QUOTED  LITERAL  OR  HOLLERITH  STRING  OR  ARRAY  = 

CONTAINING  THE  INTENDED  LABEL  FOR  THE  GRAPH.  THE  = 
MAXIMUM  ALLOWABLE  LENGTH  (INCLUDING  '$'  CHARACTER)  = 
IS  32  CHARACTERS. 

XLABEL  IS  A  QUOTED  LITERAL  OR  HOLLERITH  STRING  OR  ARRAY  = 

CONTAINING  THE  INTENDED  LABEL  OF  THE  X-AXIS  OF  THE  = 
GRAPH.  IN  THIS  PROGRAM,  XLABEL  IS  ALWAYS  'ALPHA'.  = 

YLABEL  IS  A  QUOTED  LITERAL  OR  HOLLERITH  STRING  OR  ARRAY  = 

CONTAINING  THE  INTENDED  LABEL  OF  THE  Y-AXIS  OF  THE  = 
GRAPH.  IN  THIS  PROGRAM,  YLABEL  IS  ALWAYS  'BETA*. 

MINMAX  IS  A  PARAMETER  THAT  DETERMINES  WHETHER  THE  MINIMUM  = 
AND  MAXIMUM  VALUES  FOR  THE  AXES  ARE  TO  BE  ASSIGNED  = 
BY  THE  USER,  OR  WHETHER  THEY  WILL  BE  'AUTOSCALED' .  = 

CRVTTL  IS  A  QUOTED  LITERAL  OR  HOLLERITH  STRING  OR  ARRAY  = 

AND  TERMINATED  BY  A  '$'  CHARACTER  SPECIFYING  THE  = 
INTENDED  NAME  WHICH  LABELS  AN  INDIVIDUAL  CURVE. 

CRVNUM  IS  A  REAL  VARIABLE  OR  CONSTANT  THAT  SPECIFIES  THE  = 
VALUE  TO  BE  CONCATENATED  ONTO  THE  END  OF  'CRVTTL'.  = 
FOR  EXAMPLE,  IF  THIS  CURVE  REPRESENTS  'ZETA  =  0.5'  = 

THEN  CRVTTL  =  'Z=$',  WHILE  CRVNUM  =  0.5. 


SUBROUTINE  PLOTD(XDATA,  YDATA,  NNPTS,  EJECT,  LABEL,  XLABEL, 
X  YLABEL,  MINMAX,  CRVTTL,  CRVNUM, 

X  XMIN,  XMAX,  YMIN,  YMAX , GRD) 

REALX4  XDATA(l),  YDATA(l) 

REAL*4  XMIN,  XMAX,  YMIN,  YMAX 

REALX4  CRVTTL, CRVNUM 

INTEGER*4  NPTS, NNPTS 

LQGICALX4  EJECT 

LOGICALX1  LABEL! 1 ) 

LOGICALX1  XLABEL ( 1 ) 

LQGICALXl  YLABEL ( 1 ) 

C 

L0GICALX4  INIT  /.FALSE./ 


SET  THE  NUMBER  OF  POINTS 
NPTS  =  IABS(NNPTS) 

IF  THE  ROUTINE  HAS  BEEN  NOT  BEEN  INITIALIZED  (INIT  =  .FALSE.) 
THEN  INITIALIZE  IT. 


IF  (.NOT.  INIT)  CALL  PLDOOKXDATA,  YDATA,  NPTS,  LABEL,  XLABEL, 
X  YL ABEL ,  MINMAX,  CRVTTL , CRVNUM, 

X  XMIN,  XMAX,  YMIN,  YMAX,GRD) 


FRAME  THE  PLOT  AND  REORDER  SYMBOLS 


IF  ( .NOT.  INIT)  CALL  FRAME 


INDICATE  INITIALIZATION  IN  CASE  THIS  ROUTINE  IS  RE-ENTERED 
WITH  MULTIPLE  CURVES. 


INIT  =  .TRUE. 


CALCULATE  THE  NUMBER  OF  POINTS  TO  GIVE  ABOUT  5  MARKERS  PER 
LINE 


TENTATIVELY  SET  POINTS  ONLY,  THEN  CHECK 
NMARK  =  -1 

IF( NNPTS . L E . 0 )GO  TO  10 


CURVE  WANTED,  SET  MARKERS 
IF  ( MOD( NPTS ,  4)  . EQ .  1)  NMARK  =  NPTS  /  4 
IF  ( MOD( NPTS ,  4)  .NE.  1)  NMARK  =  NPTS  /  3 
IF  (NMARK  .EQ.  0)  NMARK  =  1 
10  CONTINUE 


DRAW  THE  CURVE 

XXMAX  =  -1.0E75 
YYMAX  =  -1.0E75 

DO  20  1=1, NPTS 

IF(XDATA(I)  .GT.  XMAX  .OR.  XDATA(I)  .LT.  XMIN  .OR. 

X  YDATA(I)  .GT.  YMAX  .OR.  YDATA(I)  .LT.  YMIN)  GO  TO  20 

I F( XXMAX  .LT.  XDATA(I))  J=I 
I F( XXMAX  .LT.  XDATA(I))  XXMAX  =  XDATA(I) 

I F( YYMAX  .LT.  YDATA(I))  K=I 
I F( YYMAX  .LT.  YDATA(I))  YYMAX  =  YDATA(I) 

20  CONTINUE 

I F ( ( YMAX-YYMAX )  -  (XMAX-XXMAX) )  40,30,30 
30  L  =  J 

GO  TO  50 
40  L  =K 

50  CALL  GRACE( 0 . ) 

CALL  DOT 

IF  (GRD)  60,60,70 
60  CALL  GRID( 1,1) 

70  CALL  RESET ( 'DASH' ) 

IF  (CRVNUM  .EQ.  -9.7531)  GO  TO  80 
CALL  HEIGHT (0.125) 

CALL  RLMESS(CRVTTL , 3, XDATA( L ) ,YDATA(L  ) ) 

CALL  RL REAL (CRVNUM, 2, ' ABUT ABUT ' ) 

CALL  THKCRV(0. 015) 

80  CALL  CURVE(XDATA,  YDATA,  NPTS,  NMARK) 

CALL  RESET ( ' THKCRV ' ) 


IF  THIS  IS  NOT  THE  LAST  (OR  ONLY)  CURVE  ON  THIS  GRAPH,  THEN 


ooo  onooooo  on  oon  nnon  non  oo  oon  oo  ooo  ooo 
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EXIT.  OTHERWISE  CLOSE  THE  PLOT  AND  TURN  OFF  INITIALIZATION 
FLAG. 

IF  (.NOT.  EJECT)  RETURN 

END  OF  THIS  PLOT 

CALL  ENDPL(O) 

INIT  =  .FALSE. 

RETURN 

END 


SUBROUTINE  PLD001 ( XDATA ,  YDATA,  NPTS,  LABEL,  XLABEL , 
X  YL ABEL ,  MINMAX,  CRVTTL ,  CRVNUM, 

X  XMIN,  XMAX,  YMIN,  YMAX ,  GRD) 

THIS  SUBROUTINE  DOES  THE  INITIALIZING. 


REALX4 

XDATA(NPTS) 

REAL*4 

YDATA(NPTS) 

INTEGERX4 

NPTS 

LOGICALXl 

LABEL ( 1 ) 

LOGICALX1 

XLABEL ( 1 ) 

LOGICALXl 

YLABEL ( 1 ) 

REALX4 

XMIN 

REALX4 

XMAX 

REALX4 

YMIN 

REALX4 

YMAX 

INITIALIZE  DISSPLA 
CALL  PLD009 

CALL  HEADINCLABEL,  100,  2.,  1) 
CALL  XNAME(XL ABEL ,  100) 

CALL  YNAMECYL ABEL ,  100) 


EXTRACT  MINIMA  AND  MAXIMA 

IF  (MINMAX  .NE.  1)  GO  TO  90 

CALL  PL D01 0 ( XDAT A ,  NPTS,  XMIN,  XMAX) 

CALL  PLD010(YDATA,  NPTS,  YMIN,  YMAX) 

CALL  THE  LINEAR-LINEAR  INITIALIZING  ROUTINE 

90  CALL  PLD011  ( XMIN, XMAX, YMIN, YMAX) 

RETURN 

END 


SUBROUTINE  PLD009 

THIS  SUBROUTINE  ESTABLISHES  THE  PARAMETERS  FOR  DISSPLA. 

NOTE  THAT  IT  IS  THE  USER’S  RESPONSIBILITY  TO  NOMINATE  THE 
GRAPHIC  DEVICE. 


CALL  NOCHEK 
CALL  NOBRDR 
CALL  PAGE( 14 . , 14 . ) 
CALL  PHYS0R( 2 . , .75) 

GO  TO  LEVEL  2. 

CALL  AREA2D( 9  .  ,  9  .  ) 
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LETTERING  IS  DUPLX  WITH  UPPER  CASE  ONLY. 

CALL  DUPLX 

CALL  BASALF( 'STANDARD* ) 

INTEGER  (OR  ROUNDED)  AXES  WITH  Y  AXIS  LABELLING  AT  0  DEGREES. 

CALL  INTAXS 
CALL  YAXANGt 0 . ) 

RETURN 

END 


SUBROUTINE  PLD010(V,  N,  MIN,  MAX) 

THIS  SUBROUTINE  SCANS  A  VECTOR  FOR  MAXIMUM  AND  MINIMUM 
INPUT  PARAMETERS! 

V  DATA  VECTOR  (REAL) 

N  NUMBER  OF  POINTS  IN  VECTOR  V  (INTEGER) 

OUTPUT  PARAMETERS: 

MIN  VECTOR  ORIGIN  (REAL) 

MAX  VECTOR  MAXIMUM  (REAL) 

REAL**  V(N) 

INTEGER**  N 
REAL**  MIN 

REAL**  MAX 

INITIALIZE  THE  MAXIMA  AND  MINIMA 

MIN  =  1.0E75 
MAX  =  -1.0E75 

FIND  MAXIMUM  AND  MINIMUM  OF  VECTOR  V 

DO  100  I  =  1,  N 
IF  (MIN  .GT.  V(I))  MIN  =  V(I) 

IF  (MAX  .LT.  V( I ) )  MAX  =  V(I) 

100  CONTINUE 
RETURN 
END 


SUBROUTINE  PLD011  (XMIN, XMAX, YMIN, YMAX) 

THIS  SUBROUTINE  SETS  UP  DISSPLA  FOR  A  LINEAR-LINEAR  AXIS  PLOT 

REAL**  XMIN 

REAL**  XMAX 

REAL**  YMIN 

REAL**  YMAX 

A  SIMPLE  CALL  TO  GRAF  WILL  DO  IT... 

CALL  HEIGHT( 0 . 175 ) 

CALL  GRAF(XMIN »  'SCALE',  XMAX,  YMIN,  'SCALE',  YMAX) 

RETURN 

END 


oooooooo  ooooooo 


SUBROUTINE  ROOTS  --  CALCULATES  ROOTS  OF  THE  NO  ORDER  EQUATION 
FOR  EVERY  TENTH  VALUE  ALPHA/BETA  PAIR. 


SUBROUTINE  ROOTS  (ALPHA,  BETA,  AJ,  BJ,  CJ,  N02) 

INTEGERX4  N02,  NN,  NNN,  NNNN 

REALX4  ALPHA,  BETA,  AJ(IOO),  BJ(IOO),  CJ(IOO) 

REALX8  COEFCIOO),  ROOTMYUOO) 

WRITE(6 , 40 ) 

NN  =  N02+1 
NNNN  =  N02 
DO  10  1=1, NN 
NNN  =  NN+l-I 

10  COEFCNNN)  =  CJ(I)  +  (ALPHA  x  AJ(I))  +  (BETA  X  BJ(I)) 
CALL  ZPOLR  (COEF, N02, ROOTMY ,  IER) 

WRITE  (6,50) 

DO  20  1=1, NNNN 

II  =  2X1 

III  =  II-l 

20  WRITE  (6,60)  ROOTMY(III),  ROOTMY(II) 

WRITE  (6,30) 

30  FORMAT (////////////////////  ) 

40  FORMAT (/20X, 1  ROOTS  FOR  ABOVE  ALPHA,  BETA') 

50  F0RMAT(21X, '  REAL  PART  IMAG  PART') 

60  FORMAT (21X, E10.4,6X,E10.4) 

RETURN 

END 


SUBROUTINE  SAVIT  —  SAVES  DATA  IN  FN  FT  FM  =  INAME  DATA  Al, 
WHERE  INAME  IS  THE  USER'S  CHOICE. 


SUBROUTINE  SAVIT 

COMMON  /SAVE/  LABEL,  WN,  ND,  N02,  NC,  CJ ,  AJ,  BJ, 

X  N2,  Z ETA,  NS,  SIGMA,  NW,  W,  NZW,  ZW, 

X  XMIN,  XMAX,  YMIN,  YMAX 

REALX4  WN,  CJ(IOO),  AJ(IOO),  BJ(IOO),  XMIN,  XMAX,  YMIN,  YMAX 
REALX4  ZETA(IOO),  SIGMA(IOO),  W(100),  ZW(IOO) 

INTEGER  ND,  N02,  NC,  NZ,  NS,  NW,  NZW 
CHARACTERX4  LABEL(9),  INAME(2) 

WRITE(6 , 10 ) 


READ( 5, 20 )  (INAME(I),  1=1,2) 
CALL  FRTCMS( ' FILEDEF  ','02 

', 'DISK 

',  INAME, 'DATA 

• ) 

WRITE( 2, 30 )  ( LABEL ( I ) ,  1=1, 
WRITE(2,X)  WN 

WRITE(2, *)  ND,  N02,  NC,  NZ, 

9) 

NS, 

NW,  NZW 

WRITE(2,X)  (CJ ( J ) ,  J=l,  NC) 

WRITE(2,X)  (AJ(J) ,  J=l,  NC) 

WRITE(2,X)  ( BJ ( J ) ,  J=l,  NC) 

WRITE(2, X)  (ZETA(M) ,  M=l,  NZ) 

WRITE(2,x)  (SIGMA(M) ,  M=l,  NS) 

WRITE(2,X)  (W(M) ,  M=1 ,  NW) 

WRITE( 2 , X )  (ZW(M) ,  M=1 ,  NZW) 

WRITE(2, x)  XMIN,  XMAX,  YMIN,  YMAX 
10  F0RMAT(5X, 'UNDER  WHAT  NAME  DO  YOU  WANT  TO  SAVE  THE  DATA?',/, 
X  5X, ' (8  CHARACTERS  MAX)') 

20  FORMAT ( 2A4 ) 

30  FORMAT ( 9A4 ) 

END  FILE  02 
REWIND  02 
RETURN 
END 


C  SUBROUTINE  GETIT  —  RETRIEVES  DATA  FROM  FN  FT  FM  =  INAME  DATA  Al, 
C  WHERE  INAME  IS  THE  USER'S  CHOICE. 


SUBROUTINE  GETIT 

COMMON  /SAVE/  LABEL,  WN,  ND,  N02,  NC,  CJ,  AJ,  BJ, 

X  NZ,  ZETA,  NS,  SIGMA,  NW,  W,  NZW,  ZW, 

X  XMIN,  XMAX,  YMIN,  YMAX 

REAL*4  WN,  CJC100),  AJC100),  BJC100),  XMIN,  XMAX,  YMIN,  YMAX 
REAL**  ZETAC100),  SIGMAC100),  WC100),  ZW(IOO) 

INTEGER  ND,  N02,  NC,  NZ,  NS,  NW,  NZW 
CHARACTER**  LABEL (9) 

CHARACTERX8  INAME 
CHARACTERX21  NAME 
WRITEC 6,40) 

READC  5, 50 )  INAME 

NAME  =  'STATE  '//INAME//'  DATA  *’ 

CALL  EXCMSCNAME,RC) 

IF  CRC  .EQ.  0)  GOTO  20 
WRITE(6 , 30 ) 

GOTO  10 

CALL  FRTCMSC 'FILEDEF  ','02  »,'DISK  ’, INAME, • DATA  • 

READ(2,60)  (LABEL(I),  1=1,  9) 

READC2,*)  WN 

READ(2, *)  ND,  N02,  NC,  NZ,  NS,  NW,  NZW 
READ(2,X)  (CJ(J),  J=l,  NC) 

READ(2,X)  (AJ(J),  J=l,  NC) 

READ(2,X)  (BJCJ),  J=l,  NC) 

READ(2,X)  (ZETA(M) ,  M=l,  NZ) 

READC2,*)  (SIGMA(M),  M=l,  NS) 

READC 2, X)  (WCM) ,  M=l,  NW) 

READC2, X)  (ZW(M),  M=l,  NZW) 

READC2, x)  XMIN,  XMAX,  YMIN,  YMAX 

F0RMATC5X, 'DATA  FILE  NOT  FOUND  -  TRY  ANOTHER') 

F0RMATC5X, 'UNDER  WHAT  NAME  IS  THE  DATAFILE  SAVED?  ',/, 

X  5X, ' (8  CHARACTERS  MAX)') 

F0RMATC1A8) 

F0RMATC9A4) 

END  FILE  02 
REWIND  02 
RETURN 
END 


C  SUBROUTINE  ASTER  PLACES  A  $  AT  THE  END  OF  A  CHARACTER  STRING 


SUBROUTINE  ASTERC LLINES , LINE) 

CHARACTER**  DOLLAR/ '$  '/,  LLINESC8),  LINEC9) 

DO  10  1=1,8 

LINE(I)=LLINESCI) 

CONTINUE 

LINE(9)=D0LLAR 

RETURN 


'■  «  *  V 


SUBROUTINE  READR  --  INTERACTIVELY  READS  A  REAL  NUMBER  REPLY. 
IF  THE  USER  INADVERTENTLY  ENTERS  A  NULL  STRING 
A  WARNING  IS  ISSUED  AND  ONE  RECOVERY  IS  ALLOWED. 


SUBROUTINE  READR  (ANSR) 

REALX4  ANSR 

INTEGER  COUNT 

COUNT=0 

CONTINUE 

COUNT =COUNT+l 

IF  (COUNT.LT. 3)  GO  TO  20 

WRITE  (5,60) 

GO  TO  40 
CONTINUE 

READ  (5,x,END=30,ERR=30)  ANSR 

RETURN 

REWIND  5 

WRITE  (5,50) 

GO  TO  10 
CONTINUE 
STOP 

FORMAT  (IX,*  WARNING.  NULL  STRINGS  ARE  NOT  ALLOWED,  ENTER  A  NUMER 
XICAL  VALUE. * ) 

FORMAT  (///, 5X, '  PROGRAM  TERMINATION  -  TWO  NULL  STRINGS  ENTERED!') 
END 


SUBROUTINE  READI  —  INTERACTIVELY  READS  AN  INTEGER  REPLY. 

IF  THE  USER  INADVERTENTLY  ENTERS  A  NULL  STRING  OR  NEGATIVE  VALUE 
A  WARNING  IS  ISSUED  AND  ONE  RECOVERY  IS  ALLOWED. 


SUBROUTINE  READI  (IANS) 

INTEGER  COUNT, IANS 

COUNT=0 

CONTINUE 

C0UNT=C0UNT+1 

IF  (COUNT.LT. 3)  GO  TO  20 

WRITE  (5,70) 

GO  TO  50 


%' 

20 

CONTINUE 

READ  (5,*,END=40,ERR=40)  IANS 

IF  (IANS)  40,30,30 

1  X 
> 

30 

CONTINUE 

» 

RETURN 

L 

40 

REWIND  5 

P 

WRITE  (5,60) 

V 

L' 

GO  TO  10 

• 

V' 

50 

CONTINUE 

V 

STOP 

S' 

60 

FORMAT  (IX,*  WARNING.  IMPROPER  DATA  ENTRY!  ENTER  A  POSITIVE  INTEG 

K 

XER.') 

& 

70 

FORMAT  (/// , 5X, '  PROGRAM  TERMINATION  -  2  IMPROPER  DATA  ENTRIES!') 

f  w- 


ft 


SUBROUTINE  REA DC  —  INTERACTIVELY  READS  A  CHAR  STRING  REPLY. 
(•YES'  OR  'NO').  IF  THE  USER  INADVERTENTLY  ENTERS  A  NULL  STRING 
A  WARNING  IS  ISSUED  AND  ONE  RECOVERY  IS  ALLOWED. 


SUBROUTINE  REA DC  (CANS) 

INTEGER  COUNT 

CHARACTERS  CANS 

C0UNT=0 

CONTINUE 

COUNT =COUNT+l 

IF  (COUNT.LT. 3)  GO  TO  20 

WRITE  (5,60) 

GO  TO  AO 
CONTINUE 
REWIND  5 

READ  (5, 70, END5 30, ERR=30 )  CANS 

RETURN 

REWIND  5 

WRITE  (5,50) 

GO  TO  10 
CONTINUE 
STOP 

FORMAT  (IX,'  WARNING.  NULL  STRINGS  ARE  NOT  ALLOWED  •) 

FORMAT  (/// , 5X, '  PROGRAM  TERMINATION  -  TWO  NULL  STRINGS  ENTERED!') 
FORMAT  (A2) 


SUBROUTINE  READL  —  INTERACTIVELY  READS  A  STRING  OF  CHARACTERS.  = 

IF  THE  USER  INADVERTENTLY  ENTERS  A  NULL  STRING 
A  WARNING  IS  ISSUED  AND  ONE  RECOVERY  IS  ALLOWED. 

SUBROUTINE  READL(LLINES) 

INTEGER  COUNT,  I,  NIX 

CHARACTERS  BBLANK/'  LLINES(8) 

DO  10  1=1,8 

LLINES(I)  =  BBLANK 
10  CONTINUE 
COUNT=0 

20  COUNT =C0UNT+1 

IF(COUNT . LT . 3)  GO  TO  30 
WRITE(6,70) 

GO  TO  50 
30  CONTINUE 
REWIND  5 

READ( 5,80, END=40 , ERR=30 )(LLINES(J),J=1,9) 

RETURN 

40  REWIND  5 
WRITE(6 ,60) 

GO  TO  20 
50  CONTINUE 
STOP 

60  FORMAT ( IX, '  WARNING.  NULL  STRINGS  ARE  NOT  ALLOWED,  ENTER  CHARACTER 
X  VALUES. ') 

70  FORMAT  (///,5X,'  PROGRAM  TERMINATION  -  TWO  NULL  STRINGS  ENTERED!') 
80  FORMAT ( 9A4) 

END 


KK5W5 
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SUBROUTINE  META  --  BY  CALLING  COMPRS  AS  A  SUBROUTINE  HERE, 
DONEPL  HAS  SUFFICIENT  TIME  TO  FINISH;  OTHERWISE  COMPRS  IGNORED 
AND  METAFILE  NOT  SAVED 


SUBROUTINE  META 
DO  10  1=1,900000 
10  CONTINUE 
CALL  COMPRS 
RETURN 
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